Stick Insect Climate PCA

Data

Sorry for all of the packages. They grew and grew and grew, and I don’t want to refactor everything to reduce the dependencies.

library(raster)
library(data.table)
library(sf)
#RStoolbox has some dependencies like openMP that can be difficult to compile on a Mac (needed for the dependent package "caret"). If you have High Sierra OS or newer, search for instructions specific to your OS- it's a lot easier than older OS's.
library(RStoolbox)
library(leaflet)
library(plotly)
library(gdata)
library(BSDA)
library(ade4)
library(readxl)
library(janitor)
library(rnaturalearth)
library(adehabitatHS)
library(climateStability)
library(htmlwidgets)
library(tidyverse)
library(ggspatial)
library(here)
library(mapview)
library(fs)

# source locality plotting function
source(here::here("R", "plot_leaflet_function.R")) 

# source climate pca plotting function
source(here::here("R", "plot_climate_pca_function.R")) 

# source species pca plotting function
source(here::here("R", "species_pca_function.R")) 

# source function that creates a minimum convex polygon around points
source(here::here("R", "min_convex_poly.R")) 


source(here::here("R", "enfa_calc_function.R"))
source(here::here("R", "marginality_lollipop_plot.R"))
source(here::here("R", "presence_absence_raster_function.R"))
source(here::here("R", "crop_background_env_function.R"))
source(here::here("R", "enfa_hex_plot.R"))
source(here::here("R", "total_climate_pca_plot.R"))
source(here::here("R", "raster_correlation_function.R"))
source(here::here("R", "move_to_species.R"))

Create subdirectories (if needed) and establish paths.

# function to check if subfolder exists. If not, make it
create_dir <- function(out_path) {
  if (!dir.exists(out_path)) {
    dir.create(out_path)
  } else
    print("Directory already there.")
}

# path for interactive plots, maps, etc. output
interactive_path <- here::here("output", "interactive_plots")

create_dir(interactive_path)

[1] "Directory already there."

# climate data
climate_path <- 
  here::here("data",
             "climate")

create_dir(climate_path)

[1] "Directory already there."

# temperature climate data
temp_path <-
  here::here("data",
       "climate",
       "paleoclim_late_pleistocene",
       "temperature")

create_dir(temp_path)

[1] "Directory already there."

# precipitation climate data
precip_path <-
  here::here("data",
       "climate",
       "paleoclim_late_pleistocene",
       "precipitation")

create_dir(precip_path)

[1] "Directory already there."

# raster output
raster_path <- 
  here::here("output",
             "rasters")

create_dir(raster_path)

[1] "Directory already there."

# map (non-interactive) output
map_path <- here::here("output", "maps")

create_dir(map_path)

[1] "Directory already there."

# plot (non-interactive) output
plot_path <- here::here("output", "plots")

create_dir(plot_path)

[1] "Directory already there."

# spreadsheet output
spread_path <- here::here("output", "spreadsheets")

create_dir(spread_path)

[1] "Directory already there."

# R object (e.g. .RDS files) output
obj_path <- here::here("output", "objects")

create_dir(obj_path)

[1] "Directory already there."

Read in the spreadsheet and take a look at the data.

# data path
loc_path <- here::here("data", "all species New_6-14-19.xlsx")

# species list 
species_list <- c("horridus", 
                  "hookeri",
                  "annulata",
                  "ovobessus",
                  "huttoni",
                  "jucundum")

### read in spreadsheet
loc <- read_xlsx(loc_path) %>%
  janitor::clean_names() %>% 
  mutate(reproductive_mode = as.factor(reproductive_mode)) %>% 
  filter(species %in% species_list)

#get the number of individuals, and the sexuality counts per species
count_repro_mode <- loc %>% 
  group_by(genus, species, reproductive_mode) %>% 
  dplyr::count() %>% 
  mutate(genus_species = str_c(genus, species, sep = "_"),
         genus_species = str_replace_all(genus_species, " ", "_"),
         genus_species = str_replace_all(genus_species, "\\.", "")) %>% 
  ungroup() %>% 
  mutate(genus_species = fct_reorder(genus_species, n, sum)) %>% 
  ggplot(aes(x = genus_species, y = n, fill = reproductive_mode)) +
  geom_col() +
  coord_flip() + 
  theme_minimal()

count_repro_mode

Map

Plot a leaflet map of the localities. The leaflet map is interactive. You can click on the localities and a flag with some metadata will pop up!

#make locality shape file and assign WGS coord system
coord_points <- st_as_sf(loc, coords = c("longitude", "latitude"), 
                         crs = 4326, agr = "constant")

#use sourced plot_locs_leaflet script to plot localities
all_plot <- plot_locs_leaflet(loc, "reproductive_mode")

all_plot

# save the map
mapview::mapshot(
  all_plot,
  url = file.path(interactive_path, "all_species_map.html"),
  file = file.path(interactive_path, "all_species_map.pdf")
)

PCA-Genera

Climate Data

Read in the bioclim layers for analysis. I’m using all 19 for this preliminary exploration. I’m using CHELSA v1.2 data downloaded from their website. Plotting the first temperature layer to take a look at the data.

PCA by locality

This is a PCA of the climate data extracted for each locality, rather than a PCA of the total climate space. This gives us a general idea of differences in environmental niche.

Run the pca and check out variable loadings and proportion of variance explained by components.

Two plots: One plot of the PCA colored according to genus, with convex hulls surrounding the genera. It looks like this reflects a latitudinal gradient in temperature! You can interact with the PCA plot by clicking on points to view associated metadata. You can isolate the genus you want to view by double clicking the genus in the legend! You can also remove a genus by clicking on it once. There’s some other functionality you can explore in the toolbar at the top of the plot. The second plot is a PCA colored according to reproductive mode. It looks like asexual populations occupy slightly larger niche space, but both reproductive modes have a similar niche center.

Stability

Examining whether asexual populations occupy more unstable climates than sexual populations. Only using species with multiple sexual and asexual populations. Asexual pops tend to occupy more stable temperature environments, but less stable precipitation environments. We’re estimating stability using the method presented by Owens and Guralnick 2019- climateStability: AN R PACKAGE TO ESTIMATE CLIMATE STABILITY FROM TIME-SLICE CLIMATOLOGIES.

### Creating temperature and precipitation stability layers
### Using bio1 (average temp) and bio12 (average precip)
### 2.5 arcmin resolution, already cropped to NZ to speed up computation time


# time slices are calculated as the difference between the midpoints of the two time periods the climate layers were calculated for (e.g. midpoint of LH = (4.2 ka - 0.3 ka) / 2 + 0.3 ka = 2.25, midpoint of MH = (8.326 ka - 4.2 ka) / 2 + 4.2 = 6.263. time_slice = 6.263 - 2.25 = 4.013)

# c = current
# t_lh = midpoint time for late holocene
# t_mh = mid-Holocene
# t_eh = early-Holocene
# t_yds = Younger Dryas Stadial
# t_ba = Bølling-Allerød
# t_hs = Heinrich Stadial 1
# t_lgm = Last Glacial Maximum
c <- (2013 - 1979) / 2
t_lh <- ((4.2 - 0.3) / 2 + 0.3) * 1000
t_mh <- ((8.326 - 4.2) / 2 + 4.2) * 1000
t_eh <- ((11.7 - 8.326) / 2 + 8.326) * 1000
t_yds <- ((12.9 - 11.7) / 2 + 11.7) * 1000
t_ba <- ((14.7 - 12.9) / 2 + 12.9) * 1000
t_hs <- ((17 - 14.7) / 2 + 14.7) * 1000
t_lgm <- 21000

time_slices <- c(t_lh - c, 
                 t_mh - t_lh, 
                 t_eh - t_mh, 
                 t_yds - t_eh, 
                 t_ba - t_yds, 
                 t_hs - t_ba,
                 t_lgm - t_hs)


precip_deviation <-
  climateStability::deviationThroughTime(variableDirectory = precip_path, timeSlicePeriod = time_slices)

precip_stability <- (1 / precip_deviation)
precip_stability_scaled <- climateStability::rescale0to1(precip_stability)

# write precipitation stability to file
writeRaster(
  precip_stability_scaled,
  filename = file.path(raster_path, "precip_stability_scaled.tif"),
  format = "GTiff",
  overwrite = TRUE
)


temp_deviation <-
  climateStability::deviationThroughTime(variableDirectory = temp_path, timeSlicePeriod = time_slices)

temp_stability <- (1 / temp_deviation)
temp_stability_scaled <- climateStability::rescale0to1(temp_stability)

# write temperature stability to file
writeRaster(
  temp_stability_scaled,
  filename = file.path(raster_path, "temp_stability_scaled.tif"),
  format = "GTiff",
  overwrite = TRUE
)

overall_stability_scaled <- precip_stability_scaled * temp_stability_scaled

# write overall stability to file
writeRaster(
  overall_stability_scaled,
  filename = file.path(raster_path, "overall_stability_scaled.tif"),
  format = "GTiff",
  overwrite = TRUE
)

temp_stability_map <- ggplot() +
  ggspatial::layer_spatial(data = temp_stability_scaled) +
  labs(title = "Temperature stability") +
  scale_fill_viridis_c(na.value = "transparent") +
  theme_minimal()

precip_stability_map <- ggplot() +
  ggspatial::layer_spatial(data = precip_stability_scaled) +
  labs(title = "Precipitation stability") +
  scale_fill_viridis_c(na.value = "transparent") +
  theme_minimal()

overall_stability_map <- ggplot() +
  ggspatial::layer_spatial(data = overall_stability_scaled) +
  labs(title = "Overall stability") +
  scale_fill_viridis_c(na.value = "transparent") +
  theme_minimal()

temp_stability_map
precip_stability_map
overall_stability_map

ggsave("temp_stability.png",
       plot = temp_stability_map,
       device = "png",
       path = map_path,
       dpi = "retina")

ggsave("precip_stability.png",
       plot = precip_stability_map,
       device = "png",
       path = map_path,
       dpi = "retina")

ggsave("overall_stability.png",
       plot = overall_stability_map,
       device = "png",
       path = map_path,
       dpi = "retina")

Plot stability across species.

# filter for relevant species and asexual reproductive mode
asexual_locs <- loc %>%
  mutate(lat_long = str_c(latitude, longitude)) %>%
  filter(
    reproductive_mode == "asexual",
    species == "horridus" |
      species == "jucundum" |
      species == "hookeri" |
      species == "annulata" |
      species == "ovobessus" |
      species == "huttoni",
    !duplicated(lat_long)
  ) %>%
  dplyr::select(longitude, latitude)

# filter for relevant species and sexual reproductive mode
sexual_locs <- loc %>%
  mutate(lat_long = str_c(latitude, longitude)) %>%
  filter(
    reproductive_mode == "sexual",
    species == "horridus" |
      species == "jucundum" |
      species == "hookeri" |
      species == "annulata" |
      species == "ovobessus" |
      species == "huttoni",
    !duplicated(lat_long)
  ) %>%
  dplyr::select(longitude, latitude)

# extract preciptitation values and bind into a new dataframe
precip_asexual <-
  raster::extract(precip_stability_scaled, asexual_locs) %>%
  enframe(name = NULL, value = "precip_stability_scaled") %>%
  mutate(reproductive_mode = "asexual")

precip_sexual <-
  raster::extract(precip_stability_scaled, sexual_locs) %>%
  enframe(name = NULL, value = "precip_stability_scaled") %>%
  mutate(reproductive_mode = "sexual")

precip_df <- bind_rows(precip_asexual, precip_sexual)

# plot precipitation stability
precip_stability_plot <- ggplot(
  data = precip_df,
  aes(x = reproductive_mode, y = precip_stability_scaled, fill = reproductive_mode)
) +
  geom_violin(width = 0.8) +
  geom_boxplot(width = 0.1,
               color = "gray",
               fill = "transparent") +
  scale_fill_viridis_d(option = "magma") +
  theme_dark()

precip_stability_plot

# extract temperature values and bind into a new data frame
temp_asexual <-
  raster::extract(temp_stability_scaled, asexual_locs) %>%
  enframe(name = NULL, value = "temp_stability_scaled") %>%
  mutate(reproductive_mode = "asexual")

temp_sexual <-
  raster::extract(temp_stability_scaled, sexual_locs) %>%
  enframe(name = NULL, value = "temp_stability_scaled") %>%
  mutate(reproductive_mode = "sexual")

temp_df <- bind_rows(temp_asexual, temp_sexual)

# plot temperature stability
temp_stability_plot <- ggplot(data = temp_df,
       aes(x = reproductive_mode, y = temp_stability_scaled, fill = reproductive_mode)) +
  geom_violin(width = 0.8) +
  geom_boxplot(width = 0.1,
               color = "gray",
               fill = "transparent") +
  scale_fill_viridis_d(option = "magma") +
  theme_dark()

temp_stability_plot

# extract overall stability values and bind into a dataframe
overall_asexual <-
  raster::extract(overall_stability_scaled, asexual_locs) %>%
  enframe(name = NULL, value = "overall_stability_scaled") %>%
  mutate(reproductive_mode = "asexual")

overall_sexual <-
  raster::extract(overall_stability_scaled, sexual_locs) %>%
  enframe(name = NULL, value = "overall_stability_scaled") %>%
  mutate(reproductive_mode = "sexual")

overall_df <- bind_rows(overall_asexual, overall_sexual)

# plot overall stability
overall_stability_plot <- ggplot(
  data = overall_df,
  aes(x = reproductive_mode, y = overall_stability_scaled, fill = reproductive_mode)
) +
  geom_violin(width = 0.8) +
  geom_boxplot(width = 0.1,
               color = "gray",
               fill = "transparent") +
  scale_fill_viridis_d(option = "magma") +
  theme_dark()

overall_stability_plot

ggsave("temp_stability_plot.png",
       plot = temp_stability_plot,
       device = "png",
       path = plot_path,
       dpi = "retina")

ggsave("precip_stability_plot.png",
       plot = precip_stability_plot,
       device = "png",
       path = plot_path,
       dpi = "retina")

ggsave("overall_stability_plot.png",
       plot = overall_stability_plot,
       device = "png",
       path = plot_path,
       dpi = "retina")

Environmental space per species

These are PCAs of environmental space for species within genera. Each climate PCA is of localities for a single genus, colored by species. I’m doing this even for genera with one species, so it’s easy to see if certain localities group together.

I’m also conducting Enviromental Niche Factor Analysis to investigate the influence of reproductive mode on niche marginality.

Argosarchus horridus

# conduct pca
summary_list_ahor <- species_pca_fun(data = loc_clim, genus = "argosarchus", species = "horridus")

# plot
ahor_plot <-
  plot_clim_pca(summary_list_ahor$loc_clim,
                summary_list_ahor$summary_pca,
                factor = "reproductive_mode")
ahor_plot

htmlwidgets::saveWidget(ahor_plot,
                        file.path(interactive_path, "ahor_pca.html"),
                        selfcontained = TRUE)

#filter localities for the focal species
ahor_loc <- loc %>%
  filter(genus == "argosarchus", species == "horridus")

#use sourced plot_locs_leaflet script to plot localities
ahor_map <- plot_locs_leaflet(ahor_loc, "reproductive_mode")

ahor_map

mapview::mapshot(
  ahor_map,
  url = file.path(interactive_path, "ahor_map.html"),
  file = file.path(interactive_path, "ahor_map.pdf")
)

summary_list_ahor$summary_pca
loadings_ahor <- summary_list_ahor$summary_pca$rotation
knitr::kable(round(loadings_ahor[, 1:3], 3)) #Table of loading scores for the first 3 PCs. 

ENFA

Now I’m going to to environmental niche factor analysis between sexual and asexual populations within the species.

#get background env't for the species
ahor_bg_env <- bg_env_crop(ahor_loc, 
                           species = "horridus",
                           environment = w, 
                           buffer = 0.5)

#enfa for the sexual species
ahor_sexual_enfa <- enfa_calc_fun(locs = ahor_loc, 
                                  species = "horridus", 
                                  reproductive_mode = "sexual", 
                                  mask_raster = ahor_bg_env)

#enfa for the asexual species
ahor_asexual_enfa <- enfa_calc_fun(locs = ahor_loc, 
                                   species = "horridus", 
                                   reproductive_mode = "asexual", 
                                   mask_raster = ahor_bg_env)


#plot the marginality scores
ahor_marginality <- marginality_lollipop(sex_marg = ahor_sexual_enfa$m, 
                    asex_marg = ahor_asexual_enfa$m,
                    full_species_name = "Argosarchus horridus")

# write scores to csvs
readr::write_csv(
  ahor_asexual_enfa$m %>% enframe(name = NULL, value = "marginality"),
  file.path(spread_path, "ahor_asexual_marginality_score.csv")
)

readr::write_csv(
  ahor_sexual_enfa$m %>% enframe(name = NULL, value = "marginality"),
  file.path(spread_path, "ahor_sexual_marginality_score.csv")
)

ggsave("ahor_marginality_lollipop.png",
       plot = ahor_marginality,
       device = "png",
       path = plot_path,
       dpi = "retina")

A couple different ways to visualize the environmental variation. 1) Scatterplot visualizations of marginality vs axis 1 of the specialization with the labels removed (they make things indiscernable). Red = occupied e-space. Gray = Background e-space. The yellow dot indicates the mean marginality (it’s not the value that is on the lollipop plot, so don’t let that confuse you). 2) ENFA histogram visualizations of marginality and specialization axes. 3) PCA of total e-space with colors corresponding to sexual vs. asexual populations.

### 1) ENFA scatterplot
#access the relevant values for plotting
ahor_asexual_df <- ahor_asexual_enfa$li %>% 
  as_tibble() %>% 
  bind_cols(pr = ahor_asexual_enfa$pr)

readr::write_csv(ahor_asexual_df, 
                 file.path(spread_path, "ahor_asexual_marginality_df.csv"))

ahor_sexual_df <- ahor_sexual_enfa$li %>% 
  as_tibble() %>% 
  bind_cols(pr = ahor_sexual_enfa$pr)

readr::write_csv(ahor_sexual_df, 
                 file.path(spread_path, "ahor_sexual_marginality_df.csv"))

#asexual
ahor_enfa_spec_asexual <- enfa_hex_plot(ahor_asexual_df, marg = Mar, spec = Spe1, repro_mode = "Asexual")

ahor_enfa_spec_asexual

ggsave("ahor_enfa_spec_asexual.png",
       plot = ahor_enfa_spec_asexual,
       device = "png",
       path = plot_path,
       dpi = "retina")

#sexual
ahor_enfa_spec_sexual <- enfa_hex_plot(ahor_sexual_df, marg = Mar, spec = Spe1, repro_mode = "Sexual")

ahor_enfa_spec_sexual

ggsave("ahor_enfa_spec_sexual.png",
       plot = ahor_enfa_spec_sexual,
       device = "png",
       path = plot_path,
       dpi = "retina")

### 2) ENFA histogram
## asexual
hist(ahor_asexual_enfa)
title(main = "Asexual", adj = 0.7, line = -12)

# write to file
png(filename = file.path(plot_path, "ahor_asexual_enfa_hist.png"))
hist(ahor_asexual_enfa)
title(main = "Asexual", adj = 0.7, line = -12)
dev.off()

## sexual
hist(ahor_sexual_enfa)
title(main = "Sexual", adj = 0.7, line = -12)

# write to file
png(filename = file.path(plot_path, "ahor_sexual_enfa_hist.png"))
hist(ahor_sexual_enfa)
title(main = "Sexual", adj = 0.7, line = -12)
dev.off()

### 3) PCA of total e-space
ahor_total_pca <- total_climate_pca_plot(bg_env = ahor_bg_env, locs = ahor_loc, genus = "Argosarchus", species = "horridus")

ahor_total_pca

ggsave("ahor_total_pca.png",
       plot = ahor_total_pca,
       device = "png",
       path = plot_path,
       dpi = "retina")

ENFA Reduced

Trying out a repeat of the analyses with reduced environmental space. Prioritizing variables that will limit their distribution (i.e. variables that represent the extremes). After correlation analysis, we’re left with BIO6, BIO13, BIO14, BIO16

# PCA of background e-space. Resulting list is a correlation heatmap (cor_heatmap), a tibble of the rasters with correlation > the cutoff (default is 0.75), and a tibble of all pairwise correlations
ahor_pca <- raster_correlation(raster_stack = ahor_bg_env)
ahor_pca$cor_heatmap
ahor_pca$top_cor %>% knitr::kable()

ggsave("ahor_pca_corr.png",
       plot = ahor_pca$cor_heatmap,
       device = "png",
       path = plot_path,
       dpi = "retina")


# write correlation data frame to file
readr::write_csv(ahor_pca$top_cor, file.path(spread_path, "ahor_top_cor.csv"))

Repeat the analysis with the reduced data set. The background environment is 0.5 degrees, a ballpark dispersal limitation for all stick insect species in this study.

w_ahor <- raster::subset(w, c("bio6", "bio13", "bio14", "bio16"))

#get background env't for the species
ahor_bg_env_subset <- bg_env_crop(ahor_loc, 
                           species = "horridus",
                           environment = w_ahor, 
                           buffer = 0.5)

#enfa for the sexual species
ahor_sexual_enfa_subset <- enfa_calc_fun(locs = ahor_loc, 
                                  species = "horridus", 
                                  reproductive_mode = "sexual", 
                                  mask_raster = ahor_bg_env_subset)

#enfa for the asexual species
ahor_asexual_enfa_subset <- enfa_calc_fun(locs = ahor_loc, 
                                   species = "horridus", 
                                   reproductive_mode = "asexual", 
                                   mask_raster = ahor_bg_env_subset)

# write marginality scores to csv
readr::write_csv(
  ahor_asexual_enfa_subset$m %>% enframe(name = NULL, value = "marginality"),
  file.path(spread_path, "ahor_subset_asexual_marginality_score.csv")
)

readr::write_csv(
  ahor_sexual_enfa_subset$m %>% enframe(name = NULL, value = "marginality"),
  file.path(spread_path, "ahor_subset_sexual_marginality_score.csv")
)

#plot the marginality scores
ahor_subset_marginality <-
  marginality_lollipop(
    sex_marg = ahor_sexual_enfa_subset$m,
    asex_marg = ahor_asexual_enfa_subset$m,
    full_species_name = "Argosarchus horridus"
  )

ahor_subset_marginality

ggsave("ahor_subset_marginality.png",
       plot = ahor_subset_marginality,
       device = "png",
       path = plot_path,
       dpi = "retina")

Visualize with reduced data set

### 1) ENFA scatterplot
#access the relevant values for plotting
ahor_asexual_df_subset <- ahor_asexual_enfa_subset$li %>%
  as_tibble() %>%
  bind_cols(pr = ahor_asexual_enfa_subset$pr)

readr::write_csv(
  ahor_asexual_df_subset,
  file.path(spread_path, "ahor_subset_asexual_marginality_df.csv")
)

ahor_sexual_df_subset <- ahor_sexual_enfa_subset$li %>%
  as_tibble() %>%
  bind_cols(pr = ahor_sexual_enfa_subset$pr)

readr::write_csv(
  ahor_sexual_df_subset,
  file.path(spread_path, "ahor_subset_sexual_marginality_df.csv")
)

#asexual. Jesus these variable names are getting long
ahor_subset_enfa_spec_asexual <-
  enfa_hex_plot(
    ahor_asexual_df_subset,
    marg = Mar,
    spec = Spe1,
    repro_mode = "Asexual"
  )

ahor_subset_enfa_spec_asexual

ggsave(
  "ahor_subset_enfa_spec_asexual.png",
  plot = ahor_subset_enfa_spec_asexual,
  device = "png",
  path = plot_path,
  dpi = "retina"
)

# sexual
ahor_subset_enfa_spec_sexual <-
  enfa_hex_plot(
    ahor_sexual_df_subset,
    marg = Mar,
    spec = Spe1,
    repro_mode = "Sexual"
  )

ahor_subset_enfa_spec_sexual

ggsave(
  "ahor_subset_enfa_spec_sexual.png",
  plot = ahor_subset_enfa_spec_sexual,
  device = "png",
  path = plot_path,
  dpi = "retina"
)

### 2) ENFA histogram
# asexual
hist(ahor_asexual_enfa_subset)
title(main = "Asexual", adj = 0.7, line = -12)

png(filename = file.path(plot_path, "ahor_subset_asexual_enfa_hist.png"))
hist(ahor_asexual_enfa_subset)
title(main = "Asexual", adj = 0.7, line = -12)
dev.off()

# sexual
hist(ahor_sexual_enfa_subset)
title(main = "Sexual", adj = 0.7, line = -12)

png(filename = file.path(plot_path, "ahor_subset_sexual_enfa_hist.png"))
hist(ahor_sexual_enfa_subset)
title(main = "Sexual", adj = 0.7, line = -12)
dev.off()

### 3) PCA of total e-space
ahor_subset_total_pca <-
  total_climate_pca_plot(
    bg_env = ahor_bg_env_subset,
    locs = ahor_loc,
    genus = "Argosarchus",
    species = "horridus"
  )

ahor_subset_total_pca

ggsave(
  "ahor_subset_total_pca.png",
  plot = ahor_subset_total_pca,
  device = "png",
  path = plot_path,
  dpi = "retina"
)

# save enfa objects and remove them from the environment. They take up a lot of memory.
saveRDS(ahor_sexual_enfa, file = file.path(obj_path, "ahor_sexual_enfa.RDS"))
rm(ahor_sexual_enfa)

saveRDS(ahor_asexual_enfa, file = file.path(obj_path, "ahor_asexual_enfa.RDS"))
rm(ahor_asexual_enfa)

saveRDS(ahor_sexual_enfa_subset,
        file = file.path(obj_path, "ahor_subset_sexual_enfa.RDS"))
rm(ahor_sexual_enfa_subset)

saveRDS(ahor_asexual_enfa_subset,
        file = file.path(obj_path, "ahor_subset_asexual_enfa.RDS"))
rm(ahor_asexual_enfa_subset)

Stability

We’re also interested in seeing if asexual populations live in more unstable climatic areas relative to sexual populations.

ahor_locs_asexual <- ahor_loc %>% 
  mutate(lat_long = str_c(latitude, longitude, sep = "_")) %>% 
  filter(reproductive_mode == "asexual",
         !duplicated(lat_long)) %>% 
  dplyr::select(longitude, latitude)

ahor_locs_sexual <- ahor_loc %>% 
  mutate(lat_long = str_c(latitude, longitude, sep = "_")) %>% 
  filter(reproductive_mode == "sexual",
         !duplicated(lat_long)) %>% 
  dplyr::select(longitude, latitude)

precip_asexual_ahor <- raster::extract(precip_stability_scaled, ahor_locs_asexual) %>% 
  enframe(name = NULL, value = "precip_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

precip_sexual_ahor <- raster::extract(precip_stability_scaled, ahor_locs_sexual) %>% 
  enframe(name = NULL, value = "precip_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

precip_df_ahor <- bind_rows(precip_asexual_ahor, precip_sexual_ahor)

readr::write_csv(precip_df_ahor, 
                 file.path(spread_path, "ahor_precip_stability_df.csv"))

ahor_precip_stability_plot <- ggplot(data = precip_df_ahor, aes(x = reproductive_mode, y = precip_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

ahor_precip_stability_plot

ggsave(
  "ahor_precip_stability.png",
  plot = ahor_precip_stability_plot,
  device = "png",
  path = plot_path,
  dpi = "retina"
)

temp_asexual_ahor <- raster::extract(temp_stability_scaled, ahor_locs_asexual) %>% 
  enframe(name = NULL, value = "temp_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

temp_sexual_ahor <- raster::extract(temp_stability_scaled, ahor_locs_sexual) %>% 
  enframe(name = NULL, value = "temp_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

temp_df_ahor <- bind_rows(temp_asexual_ahor, temp_sexual_ahor)

readr::write_csv(temp_df_ahor, 
                 file.path(spread_path, "ahor_temp_stability_df.csv"))

ahor_temp_stability_plot <- ggplot(data = temp_df_ahor, aes(x = reproductive_mode, y = temp_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

ahor_temp_stability_plot

ggsave(
  "ahor_temp_stability.png",
  plot = ahor_temp_stability_plot,
  device = "png",
  path = plot_path,
  dpi = "retina"
)

overall_asexual_ahor <- raster::extract(overall_stability_scaled, ahor_locs_asexual) %>% 
  enframe(name = NULL, value = "overall_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

overall_sexual_ahor <- raster::extract(overall_stability_scaled, ahor_locs_sexual) %>% 
  enframe(name = NULL, value = "overall_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

overall_df_ahor <- bind_rows(overall_asexual_ahor, overall_sexual_ahor)

readr::write_csv(overall_df_ahor, 
                 file.path(spread_path, "ahor_overall_stability_df.csv"))


ahor_overall_stability_plot <- ggplot(data = overall_df_ahor, aes(x = reproductive_mode, y = overall_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

ahor_overall_stability_plot

ggsave(
  "ahor_overall_stability.png",
  plot = ahor_overall_stability_plot,
  device = "png",
  path = plot_path,
  dpi = "retina"
)

Seasonality

seas_precip_asexual_ahor <- raster::extract(w$bio15, ahor_locs_asexual) %>% 
  enframe(name = NULL, value = "precip_seasonality") %>% 
  mutate(reproductive_mode = "asexual")

seas_precip_sexual_ahor <- raster::extract(w$bio15, ahor_locs_sexual) %>% 
  enframe(name = NULL, value = "precip_seasonality") %>% 
  mutate(reproductive_mode = "sexual")

seas_precip_df_ahor <- bind_rows(seas_precip_asexual_ahor, seas_precip_sexual_ahor)

readr::write_csv(seas_precip_df_ahor, 
                 file.path(spread_path, "ahor_precip_seasonality_df.csv"))

ahor_precip_seasonality_plot <- ggplot(data = seas_precip_df_ahor, aes(x = reproductive_mode, y = precip_seasonality, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

ahor_precip_seasonality_plot

ggsave(
  "ahor_precip_seasonality.png",
  plot = ahor_precip_seasonality_plot,
  device = "png",
  path = plot_path,
  dpi = "retina"
)

seas_temp_asexual_ahor <- raster::extract(w$bio4, ahor_locs_asexual) %>% 
  enframe(name = NULL, value = "temp_seasonality") %>% 
  mutate(reproductive_mode = "asexual")

seas_temp_sexual_ahor <- raster::extract(w$bio4, ahor_locs_sexual) %>% 
  enframe(name = NULL, value = "temp_seasonality") %>% 
  mutate(reproductive_mode = "sexual")

seas_temp_df_ahor <- bind_rows(seas_temp_asexual_ahor, seas_temp_sexual_ahor)

readr::write_csv(seas_temp_df_ahor, 
                 file.path(spread_path, "ahor_temp_seasonality_df.csv"))

ahor_temp_seasonality_plot <- ggplot(data = seas_temp_df_ahor, aes(x = reproductive_mode, y = temp_seasonality, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

ahor_temp_seasonality_plot

ggsave(
  "ahor_temp_seasonality.png",
  plot = ahor_temp_seasonality_plot,
  device = "png",
  path = plot_path,
  dpi = "retina"
)

Move output

Put all output into species-specific subfolders.

Asteliaphasma jucundum

#pca
summary_list_ajuc <- species_pca_fun(loc_clim, 
                                     genus ="asteliaphasma", 
                                     species = "jucundum")
#plot
ajuc_plot <- plot_clim_pca(summary_list_ajuc$loc_clim, 
                           summary_list_ajuc$summary_pca, 
                           factor = "reproductive_mode")
ajuc_plot

#if selfcontained = TRUE, you can remove the folder that gets added alongside the plot. 
htmlwidgets::saveWidget(ajuc_plot, file.path(interactive_path, "ajuc_pca.html"), selfcontained = TRUE)

#filter localities for the focal species
ajuc_loc <- loc %>% 
  filter(genus == "asteliaphasma", species == "jucundum")
  
#use sourced plot_locs_leaflet script to plot localities
ajuc_map <- plot_locs_leaflet(ajuc_loc, "reproductive_mode")

ajuc_map

# save the map
mapview::mapshot(ajuc_map, url = file.path(interactive_path, "ajuc_map.html"), file = file.path(interactive_path, "ajuc_map.pdf"))

summary_list_aste$summary_pca
loadings_aste <- summary_list_aste$summary_pca$rotation
knitr::kable(round(loadings_aste[,1:3],3)) #Table of loading scores for the first 3 PCs. 

ENFA

Now I’m going to to environmental niche factor analysis between sexual and asexual populations within the species.

#get background env't for the species
ajuc_bg_env <- bg_env_crop(aste_loc, 
                           species = "jucundum",
                           environment = w, 
                           buffer = 0.5)

#enfa for the sexual species
ajuc_sexual_enfa <- enfa_calc_fun(locs = aste_loc, 
                                  species = "jucundum", 
                                  reproductive_mode = "sexual", 
                                  mask_raster = ajuc_bg_env)

#enfa for the asexual species
ajuc_asexual_enfa <- enfa_calc_fun(locs = aste_loc, 
                                   species = "jucundum", 
                                   reproductive_mode = "asexual", 
                                   mask_raster = ajuc_bg_env)


#plot the marginality scores
ajuc_marginality <-  marginality_lollipop(sex_marg = ajuc_sexual_enfa$m, 
                    asex_marg = ajuc_asexual_enfa$m,
                    full_species_name = "Asteliaphasma jucundum")

ajuc_marginality

# write scores to csvs
readr::write_csv(
  ajuc_asexual_enfa$m %>% enframe(name = NULL, value = "marginality"),
  file.path(spread_path, "ajuc_asexual_marginality_score.csv")
)

readr::write_csv(
  ajuc_sexual_enfa$m %>% enframe(name = NULL, value = "marginality"),
  file.path(spread_path, "ajuc_sexual_marginality_score.csv")
)

ggsave("ajuc_marginality_lollipop.png",
       plot = ajuc_marginality,
       device = "png",
       path = plot_path,
       dpi = "retina")

A couple different ways to visualize the environmental variation. 1) Scatterplot visualizations of marginality vs axis 1 of the specialization with the labels removed (they make things indiscernable). Red = occupied e-space. Gray = Background e-space. 2) ENFA histogram visualizations of marginality and specialization axes. 3) PCA of total e-space with colors corresponding to sexual vs. asexual populations.

### 1) ENFA scatterplot
#access the relevant values for plotting
ajuc_asexual_df <- ajuc_asexual_enfa$li %>% 
  as_tibble() %>% 
  bind_cols(pr = ajuc_asexual_enfa$pr)

readr::write_csv(ajuc_asexual_df, 
                 file.path(spread_path, "ajuc_asexual_marginality_df.csv"))


ajuc_sexual_df <- ajuc_sexual_enfa$li %>% 
  as_tibble() %>% 
  bind_cols(pr = ajuc_sexual_enfa$pr)

readr::write_csv(ajuc_sexual_df, 
                 file.path(spread_path, "ajuc_sexual_marginality_df.csv"))

#asexual
ajuc_enfa_spec_asexual <- enfa_hex_plot(ajuc_asexual_df, marg = Mar, spec = Spe1, repro_mode = "Asexual")

ggsave("ajuc_enfa_spec_asexual.png",
       plot = ajuc_enfa_spec_asexual,
       device = "png",
       path = plot_path,
       dpi = "retina")

#sexual
ajuc_enfa_spec_sexual <- enfa_hex_plot(ajuc_sexual_df, marg = Mar, spec = Spe1, repro_mode = "Sexual")

ggsave("ajuc_enfa_spec_sexual.png",
       plot = ajuc_enfa_spec_sexual,
       device = "png",
       path = plot_path,
       dpi = "retina")

### 2) ENFA histogram
#asexual
hist(ajuc_asexual_enfa)
title(main = "Asexual", adj = 0.7, line = -12)

png(filename = file.path(plot_path, "ajuc_asexual_enfa_hist.png"))
hist(ajuc_asexual_enfa)
title(main = "Asexual", adj = 0.7, line = -12)
dev.off()

# sexual
hist(ajuc_sexual_enfa)
title(main = "Sexual", adj = 0.7, line = -12)

png(filename = file.path(plot_path, "ajuc_sexual_enfa_hist.png"))
hist(ajuc_sexual_enfa)
title(main = "Sexual", adj = 0.7, line = -12)
dev.off()

### 3) PCA of total e-space
ajuc_total_pca <- total_climate_pca_plot(bg_env = ajuc_bg_env, locs = aste_loc, genus = "Asteliophasma", species = "jucundum")

ajuc_total_pca

ggsave("ajuc_total_pca.png",
       plot = ajuc_total_pca,
       device = "png",
       path = plot_path,
       dpi = "retina")

ENFA Reduced

Trying out a repeat of the analyses with reduced environmental space. Prioritizing variables that will limit their distribution (i.e. variables that represent the extremes). After correlation analysis, we’re left with BIO5, BIO6, BIO14, BIO17.

# PCA of background e-space. Resulting list is a correlation heatmap (cor_heatmap), a tibble of the rasters with correlation > the cutoff (default is 0.75), and a tibble of all pairwise correlations
ajuc_pca <- raster_correlation(raster_stack = ajuc_bg_env)
ajuc_pca$cor_heatmap
ajuc_pca$top_cor %>% knitr::kable()

ggsave("ajuc_pca_corr.png",
       plot = ajuc_pca$cor_heatmap,
       device = "png",
       path = plot_path,
       dpi = "retina")

# write correlation data frame to file
readr::write_csv(ajuc_pca$top_cor, file.path(spread_path, "ajuc_top_cor.csv"))

Repeat the analysis with the reduced data set.

w_ajuc <- raster::subset(w, c("bio5", "bio6", "bio14", "bio17"))

#get background env't for the species
ajuc_bg_env_subset <- bg_env_crop(aste_loc, 
                           species = "jucundum",
                           environment = w_ajuc, 
                           buffer = 0.5)

#enfa for the sexual species
ajuc_sexual_enfa_subset <- enfa_calc_fun(locs = aste_loc, 
                                  species = "jucundum", 
                                  reproductive_mode = "sexual", 
                                  mask_raster = ajuc_bg_env_subset)

#enfa for the asexual species
ajuc_asexual_enfa_subset <- enfa_calc_fun(locs = aste_loc, 
                                   species = "jucundum", 
                                   reproductive_mode = "asexual", 
                                   mask_raster = ajuc_bg_env_subset)

# write marginality scores to csv
readr::write_csv(
  ajuc_asexual_enfa_subset$m %>% enframe(name = NULL, value = "marginality"),
  file.path(spread_path, "ajuc_subset_asexual_marginality_score.csv")
)

readr::write_csv(
  ajuc_sexual_enfa_subset$m %>% enframe(name = NULL, value = "marginality"),
  file.path(spread_path, "ajuc_subset_sexual_marginality_score.csv")
)

# plot the marginality scores
ajuc_subset_marginality <-
  marginality_lollipop(
    sex_marg = ajuc_sexual_enfa_subset$m,
    asex_marg = ajuc_asexual_enfa_subset$m,
    full_species_name = "Asteliaphasma jucundum"
  )

ajuc_subset_marginality

ggsave("ajuc_subset_marginality.png",
       plot = ajuc_subset_marginality,
       device = "png",
       path = plot_path,
       dpi = "retina")

Visualize with the reduced data set.

### 1) ENFA scatterplot
# access the relevant values for plotting
ajuc_asexual_df_subset <- ajuc_asexual_enfa_subset$li %>% 
  as_tibble() %>% 
  bind_cols(pr = ajuc_asexual_enfa_subset$pr)

readr::write_csv(
  ajuc_asexual_df_subset,
  file.path(spread_path, "ajuc_asexual_df_subset.csv")
)


ajuc_sexual_df_subset <- ajuc_sexual_enfa_subset$li %>% 
  as_tibble() %>% 
  bind_cols(pr = ajuc_sexual_enfa_subset$pr)

readr::write_csv(
  ajuc_sexual_df_subset,
  file.path(spread_path, "ajuc_sexual_df_subset.csv")
)

# asexual
ajuc_subset_enfa_spec_asexual <- enfa_hex_plot(ajuc_asexual_df_subset, marg = Mar, spec = Spe1, repro_mode = "Asexual")

ggsave("ajuc_subset_enfa_spec_asexual.png",
       plot = ajuc_subset_enfa_spec_asexual,
       device = "png",
       path = plot_path,
       dpi = "retina")

# sexual
ajuc_subset_enfa_spec_sexual <- enfa_hex_plot(ajuc_sexual_df_subset, marg = Mar, spec = Spe1, repro_mode = "Sexual")

ggsave("ajuc_subset_enfa_spec_sexual.png",
       plot = ajuc_subset_enfa_spec_sexual,
       device = "png",
       path = plot_path,
       dpi = "retina")

### 2) ENFA histogram
# asexual
hist(ajuc_asexual_enfa_subset)
title(main = "Asexual", adj = 0.7, line = -12)

png(filename = file.path(plot_path, "ajuc_subset_asexual_enfa_hist.png"))
hist(ajuc_asexual_enfa_subset)
title(main = "Asexual", adj = 0.7, line = -12)
dev.off()


# sexual
hist(ajuc_sexual_enfa_subset)
title(main = "Sexual", adj = 0.7, line = -12)

png(filename = file.path(plot_path, "ajuc_subset_sexual_enfa_hist.png"))
hist(ajuc_sexual_enfa_subset)
title(main = "Sexual", adj = 0.7, line = -12)
dev.off()

### 3) PCA of total e-space
ajuc_subset_total_pca <- total_climate_pca_plot(bg_env = ajuc_bg_env_subset, locs = aste_loc, genus = "Asteliaphasma", species = "jucundum")

ggsave("ajuc_subset_total_pca.png",
       plot = ajuc_subset_total_pca,
       device = "png",
       path = plot_path,
       dpi = "retina")

# save enfa objects and remove them from the environment. They take up a lot of memory.
saveRDS(ajuc_sexual_enfa, file = file.path(obj_path, "ajuc_sexual_enfa.RDS"))
rm(ajuc_sexual_enfa)

saveRDS(ajuc_asexual_enfa, file = file.path(obj_path, "ajuc_asexual_enfa.RDS"))
rm(ajuc_asexual_enfa)

saveRDS(ajuc_sexual_enfa_subset,
        file = file.path(obj_path, "ajuc_subset_sexual_enfa.RDS"))
rm(ajuc_sexual_enfa_subset)

saveRDS(ajuc_asexual_enfa_subset,
        file = file.path(obj_path, "ajuc_subset_asexual_enfa.RDS"))
rm(ajuc_asexual_enfa_subset)

Stability

We’re also interested in seeing if asexual populations live in more unstable climatic areas relative to sexual populations.

ajuc_locs_asexual <- ajuc_loc %>% 
  mutate(lat_long = str_c(latitude, longitude, sep = "_")) %>% 
  filter(reproductive_mode == "asexual",
         !duplicated(lat_long)) %>% 
  dplyr::select(longitude, latitude)

ajuc_locs_sexual <- ajuc_loc %>% 
  mutate(lat_long = str_c(latitude, longitude, sep = "_")) %>% 
  filter(reproductive_mode == "sexual",
         !duplicated(lat_long)) %>% 
  dplyr::select(longitude, latitude)

precip_asexual_ajuc <- raster::extract(precip_stability_scaled, aste_locs_asexual) %>% 
  enframe(name = NULL, value = "precip_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

precip_sexual_ajuc <- raster::extract(precip_stability_scaled, aste_locs_sexual) %>% 
  enframe(name = NULL, value = "precip_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

precip_df_ajuc <- bind_rows(precip_asexual_ajuc, precip_sexual_ajuc)

readr::write_csv(precip_df_ajuc, 
                 file.path(spread_path, "ajuc_precip_stability_df.csv"))

ajuc_precip_stability_plot <- ggplot(data = precip_df_ajuc, aes(x = reproductive_mode, y = precip_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

ajuc_precip_stability_plot

ggsave("ajuc_precip_stability.png",
       plot = ajuc_precip_stability_plot,
       device = "png",
       path = plot_path,
       dpi = "retina")

temp_asexual_ajuc <- raster::extract(temp_stability_scaled, aste_locs_asexual) %>% 
  enframe(name = NULL, value = "temp_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

temp_sexual_ajuc <- raster::extract(temp_stability_scaled, aste_locs_sexual) %>% 
  enframe(name = NULL, value = "temp_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

temp_df_ajuc <- bind_rows(temp_asexual_ajuc, temp_sexual_ajuc)

readr::write_csv(temp_df_ajuc, 
                 file.path(spread_path, "ajuc_temp_stability_df.csv"))

ajuc_temp_stability_plot <- ggplot(data = temp_df_ajuc, aes(x = reproductive_mode, y = temp_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

ajuc_temp_stability_plot

ggsave("ajuc_precip_stability.png",
       plot = ajuc_temp_stability_plot,
       device = "png",
       path = plot_path,
       dpi = "retina")

overall_asexual_ajuc <- raster::extract(overall_stability_scaled, aste_locs_asexual) %>% 
  enframe(name = NULL, value = "overall_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

overall_sexual_ajuc <- raster::extract(overall_stability_scaled, aste_locs_sexual) %>% 
  enframe(name = NULL, value = "overall_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

overall_df_ajuc <- bind_rows(overall_asexual_ajuc, overall_sexual_ajuc)

readr::write_csv(overall_df_ajuc, 
                 file.path(spread_path, "ajuc_overall_stability_df.csv"))

ajuc_overall_stability_plot <- ggplot(data = overall_df_ajuc, aes(x = reproductive_mode, y = overall_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

ajuc_overall_stability_plot

ggsave("ajuc_overall_stability.png",
       plot = ajuc_overall_stability_plot,
       device = "png",
       path = plot_path,
       dpi = "retina")

Seasonality

seas_precip_asexual_ajuc <- raster::extract(w$bio15, ajuc_locs_asexual) %>% 
  enframe(name = NULL, value = "precip_seasonality") %>% 
  mutate(reproductive_mode = "asexual")

seas_precip_sexual_ajuc <- raster::extract(w$bio15, ajuc_locs_sexual) %>% 
  enframe(name = NULL, value = "precip_seasonality") %>% 
  mutate(reproductive_mode = "sexual")

seas_precip_df_ajuc <- bind_rows(seas_precip_asexual_ajuc, seas_precip_sexual_ajuc)

readr::write_csv(seas_precip_df_ajuc, 
                 file.path(spread_path, "ajuc_precip_seasonality_df.csv"))

ajuc_precip_seasonality_plot <- ggplot(data = seas_precip_df_ajuc, aes(x = reproductive_mode, y = precip_seasonality, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

ajuc_precip_seasonality_plot

ggsave(
  "ajuc_precip_seasonality.png",
  plot = ajuc_precip_seasonality_plot,
  device = "png",
  path = plot_path,
  dpi = "retina"
)

seas_temp_asexual_ajuc <- raster::extract(w$bio4, ajuc_locs_asexual) %>% 
  enframe(name = NULL, value = "temp_seasonality") %>% 
  mutate(reproductive_mode = "asexual")

seas_temp_sexual_ajuc <- raster::extract(w$bio4, ajuc_locs_sexual) %>% 
  enframe(name = NULL, value = "temp_seasonality") %>% 
  mutate(reproductive_mode = "sexual")

seas_temp_df_ajuc <- bind_rows(seas_temp_asexual_ajuc, seas_temp_sexual_ajuc)

readr::write_csv(seas_temp_df_ajuc, 
                 file.path(spread_path, "ajuc_temp_seasonality_df.csv"))

ajuc_temp_seasonality_plot <- ggplot(data = seas_temp_df_ajuc, aes(x = reproductive_mode, y = temp_seasonality, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

ajuc_temp_seasonality_plot

ggsave(
  "ajuc_temp_seasonality.png",
  plot = ajuc_temp_seasonality_plot,
  device = "png",
  path = plot_path,
  dpi = "retina"
)

Move output

Put all output into species-specific subfolders.

Clitarchus hookeri

summary_list_choo <- species_pca_fun(loc_clim, genus = "clitarchus", species = "hookeri")
choo_plot <- plot_clim_pca(summary_list_choo$loc_clim, summary_list_choo$summary_pca, factor = "reproductive_mode")
choo_plot


htmlwidgets::saveWidget(choo_plot, file.path(interactive_path, "choo_pca.html"), selfcontained = TRUE)

#filter localities for the focal genus
choo_loc <- loc %>% 
  filter(genus == "clitarchus")
  
#use sourced plot_locs_leaflet script to plot localities
choo_map <- plot_locs_leaflet(choo_loc, "reproductive_mode")

choo_map

# save the map

mapview::mapshot(
  choo_map,
  url = file.path(interactive_path, "choo_map.html"),
  file = paste0(file.path(interactive_path, "choo_map.pdf"))
)

summary_list_choo$summary_pca
loadings_choo <- summary_list_choo$summary_pca$rotation
knitr::kable(round(loadings_choo[,1:3],3)) #Table of loading scores for the first 3 PCs. 

ENFA

Now I’m going to perform environmental niche factor analysis with sexual and asexual populations within the species.

#get background env't for the species
choo_bg_env <- bg_env_crop(
  choo_loc,
  species = "hookeri",
  environment = w,
  buffer = 0.5
)

#enfa for the sexual species
choo_sexual_enfa <- enfa_calc_fun(
  locs = choo_loc,
  species = "hookeri",
  reproductive_mode = "sexual",
  mask_raster = choo_bg_env
)

#enfa for the asexual species
choo_asexual_enfa <- enfa_calc_fun(
  locs = choo_loc,
  species = "hookeri",
  reproductive_mode = "asexual",
  mask_raster = choo_bg_env
)

# write scores to csvs
readr::write_csv(
  choo_asexual_enfa$m %>% enframe(name = NULL, value = "marginality"),
  file.path(spread_path, "choo_asexual_marginality_score.csv")
)

readr::write_csv(
  choo_sexual_enfa$m %>% enframe(name = NULL, value = "marginality"),
  file.path(spread_path, "choo_sexual_marginality_score.csv")
)

#plot the marginality scores
choo_marginality <-
  marginality_lollipop(
    sex_marg = choo_sexual_enfa$m,
    asex_marg = choo_asexual_enfa$m,
    full_species_name = "Clitarchus hookeri"
  )

choo_marginality

ggsave(
  "choo_marginality_lollipop.png",
  plot = choo_marginality,
  device = "png",
  path = plot_path,
  dpi = "retina"
)

A couple different ways to visualize the environmental variation. 1) Scatterplot visualizations of marginality vs axis 1 of the specialization with the labels removed (they make things indiscernable). Red = occupied e-space. Gray = Background e-space. 2) ENFA histogram visualizations of marginality and specialization axes. 3) PCA of total e-space with colors corresponding to sexual vs. asexual populations.

### 1) ENFA scatterplot
#access the relevant values for plotting
choo_asexual_df <- choo_asexual_enfa$li %>% 
  as_tibble() %>% 
  bind_cols(pr = choo_asexual_enfa$pr)

readr::write_csv(choo_asexual_df, 
                 file.path(spread_path, "choo_asexual_marginality_df.csv"))

choo_sexual_df <- choo_sexual_enfa$li %>% 
  as_tibble() %>% 
  bind_cols(pr = choo_sexual_enfa$pr)

readr::write_csv(choo_sexual_df, 
                 file.path(spread_path, "choo_sexual_marginality_df.csv"))

# asexual
choo_enfa_spec_asexual <- enfa_hex_plot(choo_asexual_df, marg = Mar, spec = Spe1, repro_mode = "Asexual")

ggsave("choo_enfa_spec_asexual.png",
       plot = choo_enfa_spec_asexual,
       device = "png",
       path = plot_path,
       dpi = "retina")

#sexual
choo_enfa_spec_sexual <- enfa_hex_plot(choo_sexual_df, marg = Mar, spec = Spe1, repro_mode = "Sexual")

ggsave("choo_enfa_spec_sexual.png",
       plot = choo_enfa_spec_sexual,
       device = "png",
       path = plot_path,
       dpi = "retina")

### 2) ENFA histogram
#asexual
hist(choo_asexual_enfa)
title(main = "Asexual", adj = 0.7, line = -12)

png(filename = file.path(plot_path, "choo_asexual_enfa_hist.png"))
hist(choo_asexual_enfa)
title(main = "Asexual", adj = 0.7, line = -12)
dev.off()

#sexual
hist(choo_sexual_enfa)
title(main = "Sexual", adj = 0.7, line = -12)

png(filename = file.path(plot_path, "choo_sexual_enfa_hist.png"))
hist(choo_sexual_enfa)
title(main = "Sexual", adj = 0.7, line = -12)
dev.off()

### 3) PCA of total e-space
choo_total_pca <- total_climate_pca_plot(bg_env = choo_bg_env, locs = choo_loc, genus = "Clitarchus", species = "hookeri")

choo_total_pca

ggsave("choo_total_pca.png",
       plot = choo_total_pca,
       device = "png",
       path = plot_path,
       dpi = "retina")

ENFA Reduced

Trying out a repeat of the analyses with reduced environmental space. Prioritizing variables that will limit their distribution (i.e. variables that represent the extremes). After correlation analysis, we’re left with BIO8, BIO11, BIO15, BIO17.

# PCA of background e-space. Resulting list is a correlation heatmap (cor_heatmap), a tibble of the rasters with correlation > the cutoff (default is 0.75), and a tibble of all pairwise correlations
choo_pca <- raster_correlation(raster_stack = choo_bg_env)
choo_pca$cor_heatmap
choo_pca$top_cor %>% knitr::kable()

ggsave("choo_pca_corr.png",
       plot = choo_pca$cor_heatmap,
       device = "png",
       path = plot_path,
       dpi = "retina")

# write correlation data frame to file
readr::write_csv(choo_pca$top_cor, file.path(spread_path, "choo_top_cor.csv"))

Repeat the analysis with the reduced data set.

w_choo <- raster::subset(w, c("bio8", "bio11", "bio15", "bio17"))

#get background env't for the species
choo_bg_env_subset <- bg_env_crop(choo_loc, 
                           species = "hookeri",
                           environment = w_choo, 
                           buffer = 0.5)

#enfa for the sexual species
choo_sexual_enfa_subset <- enfa_calc_fun(locs = choo_loc, 
                                  species = "hookeri", 
                                  reproductive_mode = "sexual", 
                                  mask_raster = choo_bg_env_subset)

#enfa for the asexual species
choo_asexual_enfa_subset <- enfa_calc_fun(locs = choo_loc, 
                                   species = "hookeri", 
                                   reproductive_mode = "asexual", 
                                   mask_raster = choo_bg_env_subset)

# write marginality scores to csv
readr::write_csv(
  choo_asexual_enfa_subset$m %>% enframe(name = NULL, value = "marginality"),
  file.path(spread_path, "choo_subset_asexual_marginality_score.csv")
)

readr::write_csv(
  choo_sexual_enfa_subset$m %>% enframe(name = NULL, value = "marginality"),
  file.path(spread_path, "choo_subset_sexual_marginality_score.csv")
)

# plot the marginality scores
choo_subset_marginality <-
  marginality_lollipop(
    sex_marg = choo_sexual_enfa_subset$m,
    asex_marg = choo_asexual_enfa_subset$m,
    full_species_name = "Clitarchus hookeri"
  )

choo_subset_marginality

ggsave("choo_subset_marginality.png",
       plot = choo_subset_marginality,
       device = "png",
       path = plot_path,
       dpi = "retina")

Visualize with the reduced data set.

### 1) ENFA scatterplot
# access the relevant values for plotting
choo_asexual_df_subset <- choo_asexual_enfa_subset$li %>% 
  as_tibble() %>% 
  bind_cols(pr = choo_asexual_enfa_subset$pr)

readr::write_csv(
  choo_asexual_df_subset,
  file.path(spread_path, "choo_asexual_df_subset.csv")
)


choo_sexual_df_subset <- choo_sexual_enfa_subset$li %>% 
  as_tibble() %>% 
  bind_cols(pr = choo_sexual_enfa_subset$pr)

readr::write_csv(
  choo_sexual_df_subset,
  file.path(spread_path, "choo_sexual_df_subset.csv")
)

# asexual
choo_subset_enfa_spec_asexual <- enfa_hex_plot(choo_asexual_df_subset, marg = Mar, spec = Spe1, repro_mode = "Asexual")

ggsave("choo_subset_enfa_spec_asexual.png",
       plot = choo_subset_enfa_spec_asexual,
       device = "png",
       path = plot_path,
       dpi = "retina")

# sexual
choo_subset_enfa_spec_sexual <- enfa_hex_plot(choo_sexual_df_subset, marg = Mar, spec = Spe1, repro_mode = "Sexual")

ggsave("choo_subset_enfa_spec_sexual.png",
       plot = choo_subset_enfa_spec_sexual,
       device = "png",
       path = plot_path,
       dpi = "retina")

### 2) ENFA histogram
# asexual
hist(choo_asexual_enfa_subset)
title(main = "Asexual", adj = 0.7, line = -12)

png(filename = file.path(plot_path, "choo_subset_asexual_enfa_hist.png"))
hist(choo_asexual_enfa_subset)
title(main = "Asexual", adj = 0.7, line = -12)
dev.off()


# sexual
hist(choo_sexual_enfa_subset)
title(main = "Sexual", adj = 0.7, line = -12)

png(filename = file.path(plot_path, "choo_subset_sexual_enfa_hist.png"))
hist(choo_sexual_enfa_subset)
title(main = "Sexual", adj = 0.7, line = -12)
dev.off()

### 3) PCA of total e-space
choo_subset_total_pca <- total_climate_pca_plot(bg_env = choo_bg_env_subset, locs = choo_loc, genus = "Clitarchus", species = "hookeri")

choo_subset_total_pca

ggsave("choo_subset_total_pca.png",
       plot = choo_subset_total_pca,
       device = "png",
       path = plot_path,
       dpi = "retina")

# save enfa objects and remove them from the environment. They take up a lot of memory.
saveRDS(choo_sexual_enfa, file = file.path(obj_path, "choo_sexual_enfa.RDS"))
rm(choo_sexual_enfa)

saveRDS(choo_asexual_enfa, file = file.path(obj_path, "choo_asexual_enfa.RDS"))
rm(choo_asexual_enfa)

saveRDS(choo_sexual_enfa_subset,
        file = file.path(obj_path, "choo_subset_sexual_enfa.RDS"))
rm(choo_sexual_enfa_subset)

saveRDS(choo_asexual_enfa_subset,
        file = file.path(obj_path, "choo_subset_asexual_enfa.RDS"))
rm(choo_asexual_enfa_subset)

Stability

We’re also interested in seeing if asexual populations live in more unstable climatic areas relative to sexual populations.

choo_locs_asexual <- choo_loc %>% 
  mutate(lat_long = str_c(latitude, longitude, sep = "_")) %>% 
  filter(reproductive_mode == "asexual",
         !duplicated(lat_long)) %>% 
  dplyr::select(longitude, latitude)

choo_locs_sexual <- choo_loc %>% 
  mutate(lat_long = str_c(latitude, longitude, sep = "_")) %>% 
  filter(reproductive_mode == "sexual",
         !duplicated(lat_long)) %>% 
  dplyr::select(longitude, latitude)

precip_asexual_choo <- raster::extract(precip_stability_scaled, choo_locs_asexual) %>% 
  enframe(name = NULL, value = "precip_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

precip_sexual_choo <- raster::extract(precip_stability_scaled, choo_locs_sexual) %>% 
  enframe(name = NULL, value = "precip_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

precip_df_choo <- bind_rows(precip_asexual_choo, precip_sexual_choo)

readr::write_csv(precip_df_choo, 
                 file.path(spread_path, "choo_precip_stability_df.csv"))

choo_precip_stability_plot <- ggplot(data = precip_df_choo, aes(x = reproductive_mode, y = precip_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

choo_precip_stability_plot

ggsave("choo_precip_stability.png",
       plot = choo_precip_stability_plot,
       device = "png",
       path = plot_path,
       dpi = "retina")

temp_asexual_choo <- raster::extract(temp_stability_scaled, choo_locs_asexual) %>% 
  enframe(name = NULL, value = "temp_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

temp_sexual_choo <- raster::extract(temp_stability_scaled, choo_locs_sexual) %>% 
  enframe(name = NULL, value = "temp_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

temp_df_choo <- bind_rows(temp_asexual_choo, temp_sexual_choo)

readr::write_csv(temp_df_choo, 
                 file.path(spread_path, "choo_temp_stability_df.csv"))

choo_temp_stability_plot <- ggplot(data = temp_df_choo, aes(x = reproductive_mode, y = temp_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

choo_temp_stability_plot

ggsave("choo_precip_stability.png",
       plot = choo_temp_stability_plot,
       device = "png",
       path = plot_path,
       dpi = "retina")

overall_asexual_choo <- raster::extract(overall_stability_scaled, choo_locs_asexual) %>% 
  enframe(name = NULL, value = "overall_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

overall_sexual_choo <- raster::extract(overall_stability_scaled, choo_locs_sexual) %>% 
  enframe(name = NULL, value = "overall_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

overall_df_choo <- bind_rows(overall_asexual_choo, overall_sexual_choo)

readr::write_csv(overall_df_choo, 
                 file.path(spread_path, "choo_overall_stability_df.csv"))

choo_overall_stability_plot <- ggplot(data = overall_df_choo, aes(x = reproductive_mode, y = overall_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

choo_overall_stability_plot

ggsave("choo_overall_stability.png",
       plot = choo_overall_stability_plot,
       device = "png",
       path = plot_path,
       dpi = "retina")

Seasonality

seas_precip_asexual_choo <- raster::extract(w$bio15, choo_locs_asexual) %>% 
  enframe(name = NULL, value = "precip_seasonality") %>% 
  mutate(reproductive_mode = "asexual")

seas_precip_sexual_choo <- raster::extract(w$bio15, choo_locs_sexual) %>% 
  enframe(name = NULL, value = "precip_seasonality") %>% 
  mutate(reproductive_mode = "sexual")

seas_precip_df_choo <- bind_rows(seas_precip_asexual_choo, seas_precip_sexual_choo)

readr::write_csv(seas_precip_df_choo, 
                 file.path(spread_path, "choo_precip_seasonality_df.csv"))

choo_precip_seasonality_plot <- ggplot(data = seas_precip_df_choo, aes(x = reproductive_mode, y = precip_seasonality, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

choo_precip_seasonality_plot

ggsave(
  "choo_precip_seasonality.png",
  plot = choo_precip_seasonality_plot,
  device = "png",
  path = plot_path,
  dpi = "retina"
)

seas_temp_asexual_choo <- raster::extract(w$bio4, choo_locs_asexual) %>% 
  enframe(name = NULL, value = "temp_seasonality") %>% 
  mutate(reproductive_mode = "asexual")

seas_temp_sexual_choo <- raster::extract(w$bio4, choo_locs_sexual) %>% 
  enframe(name = NULL, value = "temp_seasonality") %>% 
  mutate(reproductive_mode = "sexual")

seas_temp_df_choo <- bind_rows(seas_temp_asexual_choo, seas_temp_sexual_choo)

readr::write_csv(seas_temp_df_choo, 
                 file.path(spread_path, "choo_temp_seasonality_df.csv"))

choo_temp_seasonality_plot <- ggplot(data = seas_temp_df_choo, aes(x = reproductive_mode, y = temp_seasonality, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

choo_temp_seasonality_plot

ggsave(
  "choo_temp_seasonality.png",
  plot = choo_temp_seasonality_plot,
  device = "png",
  path = plot_path,
  dpi = "retina"
)

Move output

Put all output into species-specific subfolders.

Niveaphasma annulata

summary_list_nann <- species_pca_fun(loc_clim, genus = "niveaphasma",
                                     species = "annulata")
nann_plot <- plot_clim_pca(summary_list_nann$loc_clim, summary_list_nann$summary_pca, factor = "reproductive_mode")

nann_plot

# save the plot 
htmlwidgets::saveWidget(nann_plot, file.path(interactive_path, "nann_pca.html"), selfcontained = TRUE)

#filter localities for the focal genus
nann_loc <- loc %>% 
  filter(genus == "niveaphasma")
  
#use sourced plot_locs_leaflet script to plot localities
nann_map <- plot_locs_leaflet(nann_loc, "reproductive_mode")

Assuming "longitude" and "latitude" are longitude and latitude, respectively

nann_map

# save the map
mapview::mapshot(nann_map, url = file.path(interactive_path, "nann_map.html"), file = file.path(interactive_path, "nann_map.pdf"))

summary_list_nann$summary_pca

Importance of components:
                         PC1    PC2    PC3     PC4     PC5     PC6     PC7     PC8     PC9    PC10
Standard deviation     2.991 2.4512 1.5475 0.94020 0.73312 0.39851 0.21329 0.11142 0.06490 0.04793
Proportion of Variance 0.471 0.3162 0.1260 0.04652 0.02829 0.00836 0.00239 0.00065 0.00022 0.00012
Cumulative Proportion  0.471 0.7872 0.9133 0.95977 0.98806 0.99642 0.99882 0.99947 0.99969 0.99981
                          PC11    PC12    PC13    PC14    PC15    PC16    PC17     PC18     PC19
Standard deviation     0.03918 0.02539 0.02081 0.01849 0.01623 0.01214 0.01056 0.008283 0.006705
Proportion of Variance 0.00008 0.00003 0.00002 0.00002 0.00001 0.00001 0.00001 0.000000 0.000000
Cumulative Proportion  0.99989 0.99993 0.99995 0.99997 0.99998 0.99999 0.99999 1.000000 1.000000

loadings_nann <- summary_list_nann$summary_pca$rotation
knitr::kable(round(loadings_nann[,1:3],3)) #Table of loading scores for the first 3 PCs. 

	PC1	PC2	PC3
bio1	0.091	0.373	0.184
bio10	-0.001	0.341	0.347
bio11	0.156	0.358	0.054
bio12	0.297	-0.155	0.158
bio13	0.294	-0.152	0.172
bio14	0.297	-0.158	0.134
bio15	-0.196	0.021	0.108
bio16	0.296	-0.144	0.177
bio17	0.298	-0.159	0.145
bio18	0.300	-0.144	0.131
bio19	0.313	-0.119	0.116
bio2	-0.234	-0.233	0.273
bio3	-0.215	-0.206	0.340
bio4	-0.241	-0.238	0.232
bio5	-0.119	0.211	0.496
bio6	0.177	0.345	-0.011
bio7	-0.234	-0.242	0.251
bio8	-0.171	0.255	0.194
bio9	0.116	0.150	0.291

ENFA

Now I’m going to perform environmental niche factor analysis with sexual and asexual populations within the species.

# get background env't for the species
nann_bg_env <- bg_env_crop(
  nann_loc,
  species = "annulata",
  environment = w,
  buffer = 0.5
)

#enfa for the sexual species
nann_sexual_enfa <- enfa_calc_fun(
  locs = nann_loc,
  species = "annulata",
  reproductive_mode = "sexual",
  mask_raster = nann_bg_env
)

#enfa for the asexual species
nann_asexual_enfa <- enfa_calc_fun(
  locs = nann_loc,
  species = "annulata",
  reproductive_mode = "asexual",
  mask_raster = nann_bg_env
)

# write scores to csvs
readr::write_csv(
  nann_asexual_enfa$m %>% enframe(name = NULL, value = "marginality"),
  file.path(spread_path, "nann_asexual_marginality_score.csv")
)

readr::write_csv(
  nann_sexual_enfa$m %>% enframe(name = NULL, value = "marginality"),
  file.path(spread_path, "nann_sexual_marginality_score.csv")
)

#plot the marginality scores
nann_marginality <-
  marginality_lollipop(
    sex_marg = nann_sexual_enfa$m,
    asex_marg = nann_asexual_enfa$m,
    full_species_name = "Niveaphasma annulata"
  )

nann_marginality

ggsave(
  "nann_marginality_lollipop.png",
  plot = nann_marginality,
  device = "png",
  path = plot_path,
  dpi = "retina"
)

A couple different ways to visualize the environmental variation. 1) Scatterplot visualizations of marginality vs axis 1 of the specialization with the labels removed (they make things indiscernable). Red = occupied e-space. Gray = Background e-space. 2) ENFA histogram visualizations of marginality and specialization axes. 3) PCA of total e-space with colors corresponding to sexual vs. asexual populations.

### 1) ENFA scatterplot
#access the relevant values for plotting
nann_asexual_df <- nann_asexual_enfa$li %>% 
  as_tibble() %>% 
  bind_cols(pr = nann_asexual_enfa$pr)

readr::write_csv(nann_asexual_df, 
                 file.path(spread_path, "nann_asexual_marginality_df.csv"))

nann_sexual_df <- nann_sexual_enfa$li %>% 
  as_tibble() %>% 
  bind_cols(pr = nann_sexual_enfa$pr)

readr::write_csv(nann_sexual_df, 
                 file.path(spread_path, "nann_sexual_marginality_df.csv"))

# asexual
nann_enfa_spec_asexual <- enfa_hex_plot(nann_asexual_df, marg = Mar, spec = Spe1, repro_mode = "Asexual")

ggsave("nann_enfa_spec_asexual.png",
       plot = nann_enfa_spec_asexual,
       device = "png",
       path = plot_path,
       dpi = "retina")

#sexual
nann_enfa_spec_sexual <- enfa_hex_plot(nann_sexual_df, marg = Mar, spec = Spe1, repro_mode = "Sexual")

ggsave("nann_enfa_spec_sexual.png",
       plot = nann_enfa_spec_sexual,
       device = "png",
       path = plot_path,
       dpi = "retina")

### 2) ENFA histogram
#asexual
hist(nann_asexual_enfa)
title(main = "Asexual", adj = 0.7, line = -12)

png(filename = file.path(plot_path, "nann_asexual_enfa_hist.png"))
hist(nann_asexual_enfa)

title(main = "Asexual", adj = 0.7, line = -12)
dev.off()

quartz_off_screen 
                2 

#sexual
hist(nann_sexual_enfa)
title(main = "Sexual", adj = 0.7, line = -12)

png(filename = file.path(plot_path, "nann_sexual_enfa_hist.png"))
hist(nann_sexual_enfa)

title(main = "Sexual", adj = 0.7, line = -12)
dev.off()

quartz_off_screen 
                2 

### 3) PCA of total e-space
nann_total_pca <- total_climate_pca_plot(bg_env = nann_bg_env, locs = nann_loc, genus = "Niveaphasma", species = "annulata")

nann_total_pca

ggsave("nann_total_pca.png",
       plot = nann_total_pca,
       device = "png",
       path = plot_path,
       dpi = "retina")

ENFA Reduced

Trying out a repeat of the analyses with reduced environmental space. Prioritizing variables that will limit their distribution (i.e. variables that represent the extremes). After correlation analysis, we’re left with BIO4, BIO8, BIO9, BIO11, BIO15, BIO17

# PCA of background e-space. Resulting list is a correlation heatmap (cor_heatmap), a tibble of the rasters with correlation > the cutoff (default is 0.75), and a tibble of all pairwise correlations
nann_pca <- raster_correlation(raster_stack = nann_bg_env)
nann_pca$cor_heatmap
nann_pca$top_cor %>% knitr::kable()

Var1	Var2	value
bio6	bio11	0.9929211
bio5	bio11	0.7620957
bio5	bio8	0.7548208
bio1	bio11	0.9758130
bio1	bio6	0.9440310
bio1	bio5	0.8835979
bio10	bio11	0.9027414
bio10	bio6	0.8469003
bio10	bio8	0.7688311
bio10	bio5	0.9651367
bio10	bio1	0.9744331
bio19	bio18	0.9797213
bio12	bio18	0.9937581
bio12	bio19	0.9940107
bio13	bio18	0.9925247
bio13	bio19	0.9924016
bio13	bio12	0.9989788
bio16	bio18	0.9903177
bio16	bio19	0.9930871
bio16	bio12	0.9986837
bio16	bio13	0.9993091
bio14	bio18	0.9916471
bio14	bio19	0.9922512
bio14	bio12	0.9966334
bio14	bio13	0.9943050
bio14	bio16	0.9933774
bio17	bio18	0.9872752
bio17	bio19	0.9965362
bio17	bio12	0.9972430
bio17	bio13	0.9951214
bio17	bio16	0.9952717
bio17	bio14	0.9979879
bio2	bio3	0.9295684
bio4	bio3	0.8555918
bio4	bio2	0.9812381
bio7	bio3	0.8944362
bio7	bio2	0.9940598
bio7	bio4	0.9949823

ggsave("nann_pca_corr.png",
       plot = nann_pca$cor_heatmap,
       device = "png",
       path = plot_path,
       dpi = "retina")

# write correlation data frame to file
readr::write_csv(nann_pca$top_cor, file.path(spread_path, "nann_top_cor.csv"))

Repeat the analysis with the reduced data set.

w_nann <- raster::subset(w, c("bio4", "bio8", "bio9", "bio11", "bio15", "bio17"))

#get background env't for the species
nann_bg_env_subset <- bg_env_crop(nann_loc, 
                           species = "annulata",
                           environment = w_nann, 
                           buffer = 0.5)

#enfa for the sexual species
nann_sexual_enfa_subset <- enfa_calc_fun(locs = nann_loc, 
                                  species = "annulata", 
                                  reproductive_mode = "sexual", 
                                  mask_raster = nann_bg_env_subset)

#enfa for the asexual species
nann_asexual_enfa_subset <- enfa_calc_fun(locs = nann_loc, 
                                   species = "annulata", 
                                   reproductive_mode = "asexual", 
                                   mask_raster = nann_bg_env_subset)

# write marginality scores to csv
readr::write_csv(
  nann_asexual_enfa_subset$m %>% enframe(name = NULL, value = "marginality"),
  file.path(spread_path, "nann_subset_asexual_marginality_score.csv")
)

readr::write_csv(
  nann_sexual_enfa_subset$m %>% enframe(name = NULL, value = "marginality"),
  file.path(spread_path, "nann_subset_sexual_marginality_score.csv")
)

# plot the marginality scores
nann_subset_marginality <-
  marginality_lollipop(
    sex_marg = nann_sexual_enfa_subset$m,
    asex_marg = nann_asexual_enfa_subset$m,
    full_species_name = "Niveaphasma annulata"
  )

nann_subset_marginality

ggsave("nann_subset_marginality.png",
       plot = nann_subset_marginality,
       device = "png",
       path = plot_path,
       dpi = "retina")

Visualize with the reduced data set.

### 1) ENFA scatterplot
# access the relevant values for plotting
nann_asexual_df_subset <- nann_asexual_enfa_subset$li %>% 
  as_tibble() %>% 
  bind_cols(pr = nann_asexual_enfa_subset$pr)

readr::write_csv(
  nann_asexual_df_subset,
  file.path(spread_path, "nann_asexual_df_subset.csv")
)


nann_sexual_df_subset <- nann_sexual_enfa_subset$li %>% 
  as_tibble() %>% 
  bind_cols(pr = nann_sexual_enfa_subset$pr)

readr::write_csv(
  nann_sexual_df_subset,
  file.path(spread_path, "nann_sexual_df_subset.csv")
)

# asexual
nann_subset_enfa_spec_asexual <- enfa_hex_plot(nann_asexual_df_subset, marg = Mar, spec = Spe1, repro_mode = "Asexual")

ggsave("nann_subset_enfa_spec_asexual.png",
       plot = nann_subset_enfa_spec_asexual,
       device = "png",
       path = plot_path,
       dpi = "retina")

# sexual
nann_subset_enfa_spec_sexual <- enfa_hex_plot(nann_sexual_df_subset, marg = Mar, spec = Spe1, repro_mode = "Sexual")

ggsave("nann_subset_enfa_spec_sexual.png",
       plot = nann_subset_enfa_spec_sexual,
       device = "png",
       path = plot_path,
       dpi = "retina")

### 2) ENFA histogram
# asexual
hist(nann_asexual_enfa_subset)
title(main = "Asexual", adj = 0.7, line = -12)

png(filename = file.path(plot_path, "nann_subset_asexual_enfa_hist.png"))
hist(nann_asexual_enfa_subset)

title(main = "Asexual", adj = 0.7, line = -12)
dev.off()

quartz_off_screen 
                2 

# sexual
hist(nann_sexual_enfa_subset)
title(main = "Sexual", adj = 0.7, line = -12)

png(filename = file.path(plot_path, "nann_subset_sexual_enfa_hist.png"))
hist(nann_sexual_enfa_subset)

title(main = "Sexual", adj = 0.7, line = -12)
dev.off()

quartz_off_screen 
                2 

### 3) PCA of total e-space
nann_subset_total_pca <- total_climate_pca_plot(bg_env = nann_bg_env_subset, locs = nann_loc, genus = "Niveaphasma", species = "annulata")

nann_subset_total_pca

ggsave("nann_subset_total_pca.png",
       plot = nann_subset_total_pca,
       device = "png",
       path = plot_path,
       dpi = "retina")

# save enfa objects and remove them from the environment. They take up a lot of memory.
saveRDS(nann_sexual_enfa, file = file.path(obj_path, "nann_sexual_enfa.RDS"))
rm(nann_sexual_enfa)

saveRDS(nann_asexual_enfa, file = file.path(obj_path, "nann_asexual_enfa.RDS"))
rm(nann_asexual_enfa)

saveRDS(nann_sexual_enfa_subset,
        file = file.path(obj_path, "nann_subset_sexual_enfa.RDS"))
rm(nann_sexual_enfa_subset)

saveRDS(nann_asexual_enfa_subset,
        file = file.path(obj_path, "nann_subset_asexual_enfa.RDS"))
rm(nann_asexual_enfa_subset)

Stability

We’re also interested in seeing if asexual populations live in more unstable climatic areas relative to sexual populations.

nann_locs_asexual <- nann_loc %>% 
  mutate(lat_long = str_c(latitude, longitude, sep = "_")) %>% 
  filter(reproductive_mode == "asexual",
         !duplicated(lat_long)) %>% 
  dplyr::select(longitude, latitude)

nann_locs_sexual <- nann_loc %>% 
  mutate(lat_long = str_c(latitude, longitude, sep = "_")) %>% 
  filter(reproductive_mode == "sexual",
         !duplicated(lat_long)) %>% 
  dplyr::select(longitude, latitude)

precip_asexual_nann <- raster::extract(precip_stability_scaled, nann_locs_asexual) %>% 
  enframe(name = NULL, value = "precip_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

precip_sexual_nann <- raster::extract(precip_stability_scaled, nann_locs_sexual) %>% 
  enframe(name = NULL, value = "precip_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

precip_df_nann <- bind_rows(precip_asexual_nann, precip_sexual_nann)

readr::write_csv(precip_df_nann, 
                 file.path(spread_path, "nann_precip_stability_df.csv"))

nann_precip_stability_plot <- ggplot(data = precip_df_nann, aes(x = reproductive_mode, y = precip_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

nann_precip_stability_plot

ggsave("nann_precip_stability.png",
       plot = nann_precip_stability_plot,
       device = "png",
       path = plot_path,
       dpi = "retina")

temp_asexual_nann <- raster::extract(temp_stability_scaled, nann_locs_asexual) %>% 
  enframe(name = NULL, value = "temp_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

temp_sexual_nann <- raster::extract(temp_stability_scaled, nann_locs_sexual) %>% 
  enframe(name = NULL, value = "temp_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

temp_df_nann <- bind_rows(temp_asexual_nann, temp_sexual_nann)

readr::write_csv(temp_df_nann, 
                 file.path(spread_path, "nann_temp_stability_df.csv"))

nann_temp_stability_plot <- ggplot(data = temp_df_nann, aes(x = reproductive_mode, y = temp_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

nann_temp_stability_plot

ggsave("nann_precip_stability.png",
       plot = nann_temp_stability_plot,
       device = "png",
       path = plot_path,
       dpi = "retina")

overall_asexual_nann <- raster::extract(overall_stability_scaled, nann_locs_asexual) %>% 
  enframe(name = NULL, value = "overall_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

overall_sexual_nann <- raster::extract(overall_stability_scaled, nann_locs_sexual) %>% 
  enframe(name = NULL, value = "overall_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

overall_df_nann <- bind_rows(overall_asexual_nann, overall_sexual_nann)

readr::write_csv(overall_df_nann, 
                 file.path(spread_path, "nann_overall_stability_df.csv"))

nann_overall_stability_plot <- ggplot(data = overall_df_nann, aes(x = reproductive_mode, y = overall_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

nann_overall_stability_plot

ggsave("nann_overall_stability.png",
       plot = nann_overall_stability_plot,
       device = "png",
       path = plot_path,
       dpi = "retina")

Seasonality

seas_precip_asexual_nann <- raster::extract(w$bio15, nann_locs_asexual) %>% 
  enframe(name = NULL, value = "precip_seasonality") %>% 
  mutate(reproductive_mode = "asexual")

seas_precip_sexual_nann <- raster::extract(w$bio15, nann_locs_sexual) %>% 
  enframe(name = NULL, value = "precip_seasonality") %>% 
  mutate(reproductive_mode = "sexual")

seas_precip_df_nann <- bind_rows(seas_precip_asexual_nann, seas_precip_sexual_nann)

readr::write_csv(seas_precip_df_nann, 
                 file.path(spread_path, "nann_precip_seasonality_df.csv"))

nann_precip_seasonality_plot <- ggplot(data = seas_precip_df_nann, aes(x = reproductive_mode, y = precip_seasonality, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

nann_precip_seasonality_plot

ggsave(
  "nann_precip_seasonality.png",
  plot = nann_precip_seasonality_plot,
  device = "png",
  path = plot_path,
  dpi = "retina"
)

Saving 7.29 x 4.51 in image

seas_temp_asexual_nann <- raster::extract(w$bio4, nann_locs_asexual) %>% 
  enframe(name = NULL, value = "temp_seasonality") %>% 
  mutate(reproductive_mode = "asexual")

seas_temp_sexual_nann <- raster::extract(w$bio4, nann_locs_sexual) %>% 
  enframe(name = NULL, value = "temp_seasonality") %>% 
  mutate(reproductive_mode = "sexual")

seas_temp_df_nann <- bind_rows(seas_temp_asexual_nann, seas_temp_sexual_nann)

readr::write_csv(seas_temp_df_nann, 
                 file.path(spread_path, "nann_temp_seasonality_df.csv"))

nann_temp_seasonality_plot <- ggplot(data = seas_temp_df_nann, aes(x = reproductive_mode, y = temp_seasonality, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

nann_temp_seasonality_plot

ggsave(
  "nann_temp_seasonality.png",
  plot = nann_temp_seasonality_plot,
  device = "png",
  path = plot_path,
  dpi = "retina"
)

Saving 7.29 x 4.51 in image

Move output

Put all output into species-specific subfolders.

nann_out_int_path <- file.path(interactive_path, "niveaphasma_annulata")
nann_out_plot_path <- file.path(plot_path, "niveaphasma_annulata")
nann_out_spread_path <- file.path(spread_path, "niveaphasma_annulata")
nann_out_obj_path <- file.path(obj_path, "niveaphasma_annulata")

# move interactive
move_to_species(in_path = interactive_path,
                out_path = nann_out_int_path,
                pattern = "nann")
# move plots
move_to_species(in_path = plot_path,
                out_path = nann_out_plot_path,
                pattern = "nann")
# move spreadsheets
move_to_species(in_path = spread_path,
                out_path = nann_out_spread_path,
                pattern = "nann")
# move objects
move_to_species(in_path = obj_path,
                out_path = nann_out_obj_path,
                pattern = "nann")

summary_list_tect$summary_pca
loadings_tect <- summary_list_tect$summary_pca$rotation
knitr::kable(round(loadings_tect[,1:3],3)) #Table of loading scores for the first 3 PCs. 

Tectarchus ovobessus

summary_list_tect_ovo <-
  species_pca_fun(loc_clim,
                  genus = "tectarchus",
                  species = "ovobessus")

tect_ovo_plot <-
  plot_clim_pca(summary_list_tect_ovo$loc_clim,
                summary_list_tect_ovo$summary_pca,
                factor = "reproductive_mode")

tect_ovo_plot

# save the plot
htmlwidgets::saveWidget(tect_ovo_plot,
                        file.path(interactive_path, "tect_ovo_pca.html"),
                        selfcontained = TRUE)

#filter localities for the focal species
tect_ovo_loc <- loc %>%
  filter(genus == "tectarchus", species == "ovobessus")

#use sourced plot_locs_leaflet script to plot localities
tect_ovo_map <- plot_locs_leaflet(tect_ovo_loc, "reproductive_mode")

tect_ovo_map
# save the map
mapview::mapshot(
  tect_ovo_map,
  url = file.path(interactive_path, "tect_ovo_map.html"),
  file = file.path(interactive_path, "tect_ovo_map.pdf")
)

ENFA

Now I’m going to perform environmental niche factor analysis with sexual and asexual populations within the species.

# get background env't for the species
tect_ovo_bg_env <- bg_env_crop(
  tect_loc,
  species = "ovobessus",
  environment = w,
  buffer = 0.5
)

#enfa for the sexual species
tect_ovo_sexual_enfa <- enfa_calc_fun(
  locs = tect_loc,
  species = "ovobessus",
  reproductive_mode = "sexual",
  mask_raster = tect_ovo_bg_env
)

#enfa for the asexual species
tect_ovo_asexual_enfa <- enfa_calc_fun(
  locs = tect_loc,
  species = "ovobessus",
  reproductive_mode = "asexual",
  mask_raster = tect_ovo_bg_env
)

# write scores to csvs
readr::write_csv(
  tect_ovo_asexual_enfa$m %>% enframe(name = NULL, value = "marginality"),
  file.path(spread_path, "tect_ovo_asexual_marginality_score.csv")
)

readr::write_csv(
  tect_ovo_sexual_enfa$m %>% enframe(name = NULL, value = "marginality"),
  file.path(spread_path, "tect_ovo_sexual_marginality_score.csv")
)

#plot the marginality scores
tect_ovo_marginality <-
  marginality_lollipop(
    sex_marg = tect_ovo_sexual_enfa$m,
    asex_marg = tect_ovo_asexual_enfa$m,
    full_species_name = "Tectarchus ovobessus"
  )

tect_ovo_marginality

ggsave(
  "tect_ovo_marginality_lollipop.png",
  plot = tect_ovo_marginality,
  device = "png",
  path = plot_path,
  dpi = "retina"
)

A couple different ways to visualize the environmental variation. 1) Scatterplot visualizations of marginality vs axis 1 of the specialization with the labels removed (they make things indiscernable). Red = occupied e-space. Gray = Background e-space. 2) ENFA histogram visualizations of marginality and specialization axes. 3) PCA of total e-space with colors corresponding to sexual vs. asexual populations.

### 1) ENFA scatterplot
#access the relevant values for plotting
tect_ovo_asexual_df <- tect_ovo_asexual_enfa$li %>% 
  as_tibble() %>% 
  bind_cols(pr = tect_ovo_asexual_enfa$pr)

readr::write_csv(tect_ovo_asexual_df, 
                 file.path(spread_path, "tect_ovo_asexual_marginality_df.csv"))

tect_ovo_sexual_df <- tect_ovo_sexual_enfa$li %>% 
  as_tibble() %>% 
  bind_cols(pr = tect_ovo_sexual_enfa$pr)

readr::write_csv(tect_ovo_sexual_df, 
                 file.path(spread_path, "tect_ovo_sexual_marginality_df.csv"))

# asexual
tect_ovo_enfa_spec_asexual <- enfa_hex_plot(tect_ovo_asexual_df, marg = Mar, spec = Spe1, repro_mode = "Asexual")

ggsave("tect_ovo_enfa_spec_asexual.png",
       plot = tect_ovo_enfa_spec_asexual,
       device = "png",
       path = plot_path,
       dpi = "retina")

#sexual
tect_ovo_enfa_spec_sexual <- enfa_hex_plot(tect_ovo_sexual_df, marg = Mar, spec = Spe1, repro_mode = "Sexual")

ggsave("tect_ovo_enfa_spec_sexual.png",
       plot = tect_ovo_enfa_spec_sexual,
       device = "png",
       path = plot_path,
       dpi = "retina")

### 2) ENFA histogram
#asexual
hist(tect_ovo_asexual_enfa)
title(main = "Asexual", adj = 0.7, line = -12)

png(filename = file.path(plot_path, "tect_ovo_asexual_enfa_hist.png"))
hist(tect_ovo_asexual_enfa)
title(main = "Asexual", adj = 0.7, line = -12)
dev.off()

#sexual
hist(tect_ovo_sexual_enfa)
title(main = "Sexual", adj = 0.7, line = -12)

png(filename = file.path(plot_path, "tect_ovo_sexual_enfa_hist.png"))
hist(tect_ovo_sexual_enfa)
title(main = "Sexual", adj = 0.7, line = -12)
dev.off()

### 3) PCA of total e-space
tect_ovo_total_pca <- total_climate_pca_plot(bg_env = tect_ovo_bg_env, locs = tect_loc, genus = "Tectarchus", species = "ovobessus")

tect_ovo_total_pca

ggsave("tect_ovo_total_pca.png",
       plot = tect_ovo_total_pca,
       device = "png",
       path = plot_path,
       dpi = "retina")

Trying out a repeat of the analyses with reduced environmental space. Prioritizing variables that will limit their distribution (i.e. variables that represent the extremes). After correlation analysis, we’re left with BIO4, BIO8, BIO11, BIO15, BIO17

# PCA of background e-space. Resulting list is a correlation heatmap (cor_heatmap), a tibble of the rasters with correlation > the cutoff (default is 0.75), and a tibble of all pairwise correlations
tect_ovo_pca <- raster_correlation(raster_stack = tect_ovo_bg_env)
tect_ovo_pca$cor_heatmap
tect_ovo_pca$top_cor %>% knitr::kable()

ggsave("tect_ovo_pca_corr.png",
       plot = tect_ovo_pca$cor_heatmap,
       device = "png",
       path = plot_path,
       dpi = "retina")

# write correlation data frame to file
readr::write_csv(tect_ovo_pca$top_cor, file.path(spread_path, "tect_ovo_top_cor.csv"))

ENFA Reduced

Repeat the analysis with the reduced data set.

w_tect_ovo <- raster::subset(w, c("bio4", "bio8", "bio11", "bio15", "bio17"))

#get background env't for the species
tect_ovo_bg_env_subset <- bg_env_crop(tect_loc, 
                           species = "ovobessus",
                           environment = w_tect_ovo, 
                           buffer = 0.5)

#enfa for the sexual species
tect_ovo_sexual_enfa_subset <- enfa_calc_fun(locs = tect_loc, 
                                  species = "ovobessus", 
                                  reproductive_mode = "sexual", 
                                  mask_raster = tect_ovo_bg_env_subset)

#enfa for the asexual species
tect_ovo_asexual_enfa_subset <- enfa_calc_fun(locs = tect_loc, 
                                   species = "ovobessus", 
                                   reproductive_mode = "asexual", 
                                   mask_raster = tect_ovo_bg_env_subset)

# write marginality scores to csv
readr::write_csv(
  tect_ovo_asexual_enfa_subset$m %>% enframe(name = NULL, value = "marginality"),
  file.path(spread_path, "tect_ovo_subset_asexual_marginality_score.csv")
)

readr::write_csv(
  tect_ovo_sexual_enfa_subset$m %>% enframe(name = NULL, value = "marginality"),
  file.path(spread_path, "tect_ovo_subset_sexual_marginality_score.csv")
)

# plot the marginality scores
tect_ovo_subset_marginality <-
  marginality_lollipop(
    sex_marg = tect_ovo_sexual_enfa_subset$m,
    asex_marg = tect_ovo_asexual_enfa_subset$m,
    full_species_name = "Tectarchus ovobessus"
  )

tect_ovo_subset_marginality

ggsave("tect_ovo_subset_marginality.png",
       plot = tect_ovo_subset_marginality,
       device = "png",
       path = plot_path,
       dpi = "retina")

Visualize with the reduced data set.

### 1) ENFA scatterplot
# access the relevant values for plotting
tect_ovo_asexual_df_subset <- tect_ovo_asexual_enfa_subset$li %>% 
  as_tibble() %>% 
  bind_cols(pr = tect_ovo_asexual_enfa_subset$pr)

readr::write_csv(
  tect_ovo_asexual_df_subset,
  file.path(spread_path, "tect_ovo_asexual_df_subset.csv")
)


tect_ovo_sexual_df_subset <- tect_ovo_sexual_enfa_subset$li %>% 
  as_tibble() %>% 
  bind_cols(pr = tect_ovo_sexual_enfa_subset$pr)

readr::write_csv(
  tect_ovo_sexual_df_subset,
  file.path(spread_path, "tect_ovo_sexual_df_subset.csv")
)

# asexual
tect_ovo_subset_enfa_spec_asexual <- enfa_hex_plot(tect_ovo_asexual_df_subset, marg = Mar, spec = Spe1, repro_mode = "Asexual")

ggsave("tect_ovo_subset_enfa_spec_asexual.png",
       plot = tect_ovo_subset_enfa_spec_asexual,
       device = "png",
       path = plot_path,
       dpi = "retina")

# sexual
tect_ovo_subset_enfa_spec_sexual <- enfa_hex_plot(tect_ovo_sexual_df_subset, marg = Mar, spec = Spe1, repro_mode = "Sexual")

ggsave("tect_ovo_subset_enfa_spec_sexual.png",
       plot = tect_ovo_subset_enfa_spec_sexual,
       device = "png",
       path = plot_path,
       dpi = "retina")

### 2) ENFA histogram
# asexual
hist(tect_ovo_asexual_enfa_subset)
title(main = "Asexual", adj = 0.7, line = -12)

png(filename = file.path(plot_path, "tect_ovo_subset_asexual_enfa_hist.png"))
hist(tect_ovo_asexual_enfa_subset)
title(main = "Asexual", adj = 0.7, line = -12)
dev.off()


# sexual
hist(tect_ovo_sexual_enfa_subset)
title(main = "Sexual", adj = 0.7, line = -12)

png(filename = file.path(plot_path, "tect_ovo_subset_sexual_enfa_hist.png"))
hist(tect_ovo_sexual_enfa_subset)
title(main = "Sexual", adj = 0.7, line = -12)
dev.off()

### 3) PCA of total e-space
tect_ovo_subset_total_pca <- total_climate_pca_plot(bg_env = tect_ovo_bg_env_subset, locs = tect_loc, genus = "Tectarchus", species = "ovobessus")

tect_ovo_subset_total_pca

ggsave("tect_ovo_subset_total_pca.png",
       plot = tect_ovo_subset_total_pca,
       device = "png",
       path = plot_path,
       dpi = "retina")

# save enfa objects and remove them from the environment. They take up a lot of memory.
saveRDS(tect_ovo_sexual_enfa, file = file.path(obj_path, "tect_ovo_sexual_enfa.RDS"))
rm(tect_ovo_sexual_enfa)

saveRDS(tect_ovo_asexual_enfa, file = file.path(obj_path, "tect_ovo_asexual_enfa.RDS"))
rm(tect_ovo_asexual_enfa)

saveRDS(tect_ovo_sexual_enfa_subset,
        file = file.path(obj_path, "tect_ovo_subset_sexual_enfa.RDS"))
rm(tect_ovo_sexual_enfa_subset)

saveRDS(tect_ovo_asexual_enfa_subset,
        file = file.path(obj_path, "tect_ovo_subset_asexual_enfa.RDS"))
rm(tect_ovo_asexual_enfa_subset)

Stability

We’re also interested in seeing if asexual populations live in more unstable climatic areas relative to sexual populations.

tect_ovo_locs_asexual <- tect_ovo_loc %>% 
  mutate(lat_long = str_c(latitude, longitude, sep = "_")) %>% 
  filter(reproductive_mode == "asexual",
         species == "ovobessus",
         !duplicated(lat_long)) %>% 
  dplyr::select(longitude, latitude)

tect_ovo_locs_sexual <- tect_ovo_loc %>% 
  mutate(lat_long = str_c(latitude, longitude, sep = "_")) %>% 
  filter(reproductive_mode == "sexual",
         species == "ovobessus",
         !duplicated(lat_long)) %>% 
  dplyr::select(longitude, latitude)

precip_asexual_tect_ovo <- raster::extract(precip_stability_scaled, tect_ovo_locs_asexual) %>% 
  enframe(name = NULL, value = "precip_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

precip_sexual_tect_ovo <- raster::extract(precip_stability_scaled, tect_ovo_locs_sexual) %>% 
  enframe(name = NULL, value = "precip_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

precip_df_tect_ovo <- bind_rows(precip_asexual_tect_ovo, precip_sexual_tect_ovo)

readr::write_csv(precip_df_tect_ovo, 
                 file.path(spread_path, "tect_ovo_precip_stability_df.csv"))

tect_ovo_precip_stability_plot <- ggplot(data = precip_df_tect_ovo, aes(x = reproductive_mode, y = precip_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

tect_ovo_precip_stability_plot

ggsave("tect_ovo_precip_stability.png",
       plot = tect_ovo_precip_stability_plot,
       device = "png",
       path = plot_path,
       dpi = "retina")

temp_asexual_tect_ovo <- raster::extract(temp_stability_scaled, tect_ovo_locs_asexual) %>% 
  enframe(name = NULL, value = "temp_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

temp_sexual_tect_ovo <- raster::extract(temp_stability_scaled, tect_ovo_locs_sexual) %>% 
  enframe(name = NULL, value = "temp_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

temp_df_tect_ovo <- bind_rows(temp_asexual_tect_ovo, temp_sexual_tect_ovo)

readr::write_csv(temp_df_tect_ovo, 
                 file.path(spread_path, "tect_ovo_temp_stability_df.csv"))

tect_ovo_temp_stability_plot <- ggplot(data = temp_df_tect_ovo, aes(x = reproductive_mode, y = temp_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

tect_ovo_temp_stability_plot

ggsave("tect_ovo_precip_stability.png",
       plot = tect_ovo_temp_stability_plot,
       device = "png",
       path = plot_path,
       dpi = "retina")

overall_asexual_tect_ovo <- raster::extract(overall_stability_scaled, tect_ovo_locs_asexual) %>% 
  enframe(name = NULL, value = "overall_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

overall_sexual_tect_ovo <- raster::extract(overall_stability_scaled, tect_ovo_locs_sexual) %>% 
  enframe(name = NULL, value = "overall_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

overall_df_tect_ovo <- bind_rows(overall_asexual_tect_ovo, overall_sexual_tect_ovo)

readr::write_csv(overall_df_tect_ovo, 
                 file.path(spread_path, "tect_ovo_overall_stability_df.csv"))

tect_ovo_overall_stability_plot <- ggplot(data = overall_df_tect_ovo, aes(x = reproductive_mode, y = overall_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

tect_ovo_overall_stability_plot

ggsave("tect_ovo_overall_stability.png",
       plot = tect_ovo_overall_stability_plot,
       device = "png",
       path = plot_path,
       dpi = "retina")

Seasonality

seas_precip_asexual_tect_ovo <- raster::extract(w$bio15, tect_ovo_locs_asexual) %>% 
  enframe(name = NULL, value = "precip_seasonality") %>% 
  mutate(reproductive_mode = "asexual")

seas_precip_sexual_tect_ovo <- raster::extract(w$bio15, tect_ovo_locs_sexual) %>% 
  enframe(name = NULL, value = "precip_seasonality") %>% 
  mutate(reproductive_mode = "sexual")

seas_precip_df_tect_ovo <- bind_rows(seas_precip_asexual_tect_ovo, seas_precip_sexual_tect_ovo)

readr::write_csv(seas_precip_df_tect_ovo, 
                 file.path(spread_path, "tect_ovo_precip_seasonality_df.csv"))

tect_ovo_precip_seasonality_plot <- ggplot(data = seas_precip_df_tect_ovo, aes(x = reproductive_mode, y = precip_seasonality, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

tect_ovo_precip_seasonality_plot

ggsave(
  "tect_ovo_precip_seasonality.png",
  plot = tect_ovo_precip_seasonality_plot,
  device = "png",
  path = plot_path,
  dpi = "retina"
)

seas_temp_asexual_tect_ovo <- raster::extract(w$bio4, tect_ovo_locs_asexual) %>% 
  enframe(name = NULL, value = "temp_seasonality") %>% 
  mutate(reproductive_mode = "asexual")

seas_temp_sexual_tect_ovo <- raster::extract(w$bio4, tect_ovo_locs_sexual) %>% 
  enframe(name = NULL, value = "temp_seasonality") %>% 
  mutate(reproductive_mode = "sexual")

seas_temp_df_tect_ovo <- bind_rows(seas_temp_asexual_tect_ovo, seas_temp_sexual_tect_ovo)

readr::write_csv(seas_temp_df_tect_ovo, 
                 file.path(spread_path, "tect_ovo_temp_seasonality_df.csv"))

tect_ovo_temp_seasonality_plot <- ggplot(data = seas_temp_df_tect_ovo, aes(x = reproductive_mode, y = temp_seasonality, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

tect_ovo_temp_seasonality_plot

ggsave(
  "tect_ovo_temp_seasonality.png",
  plot = tect_ovo_temp_seasonality_plot,
  device = "png",
  path = plot_path,
  dpi = "retina"
)

Move output

Put all output into species-specific subfolders.

Tectarchus huttoni

summary_list_tect_hutt <-
  species_pca_fun(loc_clim,
                  genus = "tectarchus",
                  species = "huttoni")

tect_hutt_plot <-
  plot_clim_pca(summary_list_tect_hutt$loc_clim,
                summary_list_tect_hutt$summary_pca,
                factor = "reproductive_mode")

tect_hutt_plot

# save the plot
htmlwidgets::saveWidget(tect_hutt_plot,
                        file.path(interactive_path, "tect_hutt_pca.html"),
                        selfcontained = TRUE)

#filter localities for the focal species
tect_hutt_loc <- loc %>%
  filter(genus == "tectarchus", species == "huttoni")

#use sourced plot_locs_leaflet script to plot localities
tect_hutt_map <- plot_locs_leaflet(tect_hutt_loc, "reproductive_mode")

tect_hutt_map
# save the map
mapview::mapshot(
  tect_hutt_map,
  url = file.path(interactive_path, "tect_hutt_map.html"),
  file = file.path(interactive_path, "tect_hutt_map.pdf")
)

ENFA

Now I’m going to perform environmental niche factor analysis with sexual and asexual populations within the species.

# get background env't for the species
tect_hutt_bg_env <- bg_env_crop(
  tect_loc,
  species = "huttoni",
  environment = w,
  buffer = 0.5
)

#enfa for the sexual species
tect_hutt_sexual_enfa <- enfa_calc_fun(
  locs = tect_loc,
  species = "huttoni",
  reproductive_mode = "sexual",
  mask_raster = tect_hutt_bg_env
)

#enfa for the asexual species
tect_hutt_asexual_enfa <- enfa_calc_fun(
  locs = tect_loc,
  species = "huttoni",
  reproductive_mode = "asexual",
  mask_raster = tect_hutt_bg_env
)

# write scores to csvs
readr::write_csv(
  tect_hutt_asexual_enfa$m %>% enframe(name = NULL, value = "marginality"),
  file.path(spread_path, "tect_hutt_asexual_marginality_score.csv")
)

readr::write_csv(
  tect_hutt_sexual_enfa$m %>% enframe(name = NULL, value = "marginality"),
  file.path(spread_path, "tect_hutt_sexual_marginality_score.csv")
)

#plot the marginality scores
tect_hutt_marginality <-
  marginality_lollipop(
    sex_marg = tect_hutt_sexual_enfa$m,
    asex_marg = tect_hutt_asexual_enfa$m,
    full_species_name = "Tectarchus huttoni"
  )

tect_hutt_marginality

ggsave(
  "tect_hutt_marginality_lollipop.png",
  plot = tect_hutt_marginality,
  device = "png",
  path = plot_path,
  dpi = "retina"
)

A couple different ways to visualize the environmental variation. 1) Scatterplot visualizations of marginality vs axis 1 of the specialization with the labels removed (they make things indiscernable). Red = occupied e-space. Gray = Background e-space. 2) ENFA histogram visualizations of marginality and specialization axes. 3) PCA of total e-space with colors corresponding to sexual vs. asexual populations.

### 1) ENFA scatterplot
#access the relevant values for plotting
tect_hutt_asexual_df <- tect_hutt_asexual_enfa$li %>% 
  as_tibble() %>% 
  bind_cols(pr = tect_hutt_asexual_enfa$pr)

readr::write_csv(tect_hutt_asexual_df, 
                 file.path(spread_path, "tect_hutt_asexual_marginality_df.csv"))

tect_hutt_sexual_df <- tect_hutt_sexual_enfa$li %>% 
  as_tibble() %>% 
  bind_cols(pr = tect_hutt_sexual_enfa$pr)

readr::write_csv(tect_hutt_sexual_df, 
                 file.path(spread_path, "tect_hutt_sexual_marginality_df.csv"))

# asexual
tect_hutt_enfa_spec_asexual <- enfa_hex_plot(tect_hutt_asexual_df, marg = Mar, spec = Spe1, repro_mode = "Asexual")

ggsave("tect_hutt_enfa_spec_asexual.png",
       plot = tect_hutt_enfa_spec_asexual,
       device = "png",
       path = plot_path,
       dpi = "retina")

#sexual
tect_hutt_enfa_spec_sexual <- enfa_hex_plot(tect_hutt_sexual_df, marg = Mar, spec = Spe1, repro_mode = "Sexual")

ggsave("tect_hutt_enfa_spec_sexual.png",
       plot = tect_hutt_enfa_spec_sexual,
       device = "png",
       path = plot_path,
       dpi = "retina")

### 2) ENFA histogram
#asexual
hist(tect_hutt_asexual_enfa)
title(main = "Asexual", adj = 0.7, line = -12)

png(filename = file.path(plot_path, "tect_hutt_asexual_enfa_hist.png"))
hist(tect_hutt_asexual_enfa)
title(main = "Asexual", adj = 0.7, line = -12)
dev.off()

#sexual
hist(tect_hutt_sexual_enfa)
title(main = "Sexual", adj = 0.7, line = -12)

png(filename = file.path(plot_path, "tect_hutt_sexual_enfa_hist.png"))
hist(tect_hutt_sexual_enfa)
title(main = "Sexual", adj = 0.7, line = -12)
dev.off()

### 3) PCA of total e-space
tect_hutt_total_pca <- total_climate_pca_plot(bg_env = tect_hutt_bg_env, locs = tect_loc, genus = "Tectarchus", species = "huttoni")

tect_hutt_total_pca

ggsave("tect_hutt_total_pca.png",
       plot = tect_hutt_total_pca,
       device = "png",
       path = plot_path,
       dpi = "retina")

ENFA Reduced

Trying out a repeat of the analyses with reduced environmental space. Prioritizing variables that will limit their distribution (i.e. variables that represent the extremes). After correlation analysis, we’re left with BIO4, BIO8, BIO11, BIO15, BIO17

# PCA of background e-space. Resulting list is a correlation heatmap (cor_heatmap), a tibble of the rasters with correlation > the cutoff (default is 0.75), and a tibble of all pairwise correlations
tect_hutt_pca <- raster_correlation(raster_stack = tect_hutt_bg_env)
tect_hutt_pca$cor_heatmap
tect_hutt_pca$top_cor %>% knitr::kable()

ggsave("tect_hutt_pca_corr.png",
       plot = tect_hutt_pca$cor_heatmap,
       device = "png",
       path = plot_path,
       dpi = "retina")

# write correlation data frame to file
readr::write_csv(tect_hutt_pca$top_cor, file.path(spread_path, "tect_hutt_top_cor.csv"))

Repeat the analysis with the reduced data set.

w_tect_hutt <- raster::subset(w, c("bio4", "bio8", "bio11", "bio15", "bio17"))

#get background env't for the species
tect_hutt_bg_env_subset <- bg_env_crop(tect_loc, 
                           species = "huttoni",
                           environment = w_tect_hutt, 
                           buffer = 0.5)

#enfa for the sexual species
tect_hutt_sexual_enfa_subset <- enfa_calc_fun(locs = tect_loc, 
                                  species = "huttoni", 
                                  reproductive_mode = "sexual", 
                                  mask_raster = tect_hutt_bg_env_subset)

#enfa for the asexual species
tect_hutt_asexual_enfa_subset <- enfa_calc_fun(locs = tect_loc, 
                                   species = "huttoni", 
                                   reproductive_mode = "asexual", 
                                   mask_raster = tect_hutt_bg_env_subset)

# write marginality scores to csv
readr::write_csv(
  tect_hutt_asexual_enfa_subset$m %>% enframe(name = NULL, value = "marginality"),
  file.path(spread_path, "tect_hutt_subset_asexual_marginality_score.csv")
)

readr::write_csv(
  tect_hutt_sexual_enfa_subset$m %>% enframe(name = NULL, value = "marginality"),
  file.path(spread_path, "tect_hutt_subset_sexual_marginality_score.csv")
)

# plot the marginality scores
tect_hutt_subset_marginality <-
  marginality_lollipop(
    sex_marg = tect_hutt_sexual_enfa_subset$m,
    asex_marg = tect_hutt_asexual_enfa_subset$m,
    full_species_name = "Tectarchus huttoni"
  )

tect_hutt_subset_marginality

ggsave("tect_hutt_subset_marginality.png",
       plot = tect_hutt_subset_marginality,
       device = "png",
       path = plot_path,
       dpi = "retina")

Visualize with the reduced data set.

### 1) ENFA scatterplot
# access the relevant values for plotting
tect_hutt_asexual_df_subset <- tect_hutt_asexual_enfa_subset$li %>% 
  as_tibble() %>% 
  bind_cols(pr = tect_hutt_asexual_enfa_subset$pr)

readr::write_csv(
  tect_hutt_asexual_df_subset,
  file.path(spread_path, "tect_hutt_asexual_df_subset.csv")
)


tect_hutt_sexual_df_subset <- tect_hutt_sexual_enfa_subset$li %>% 
  as_tibble() %>% 
  bind_cols(pr = tect_hutt_sexual_enfa_subset$pr)

readr::write_csv(
  tect_hutt_sexual_df_subset,
  file.path(spread_path, "tect_hutt_sexual_df_subset.csv")
)

# asexual
tect_hutt_subset_enfa_spec_asexual <- enfa_hex_plot(tect_hutt_asexual_df_subset, marg = Mar, spec = Spe1, repro_mode = "Asexual")

ggsave("tect_hutt_subset_enfa_spec_asexual.png",
       plot = tect_hutt_subset_enfa_spec_asexual,
       device = "png",
       path = plot_path,
       dpi = "retina")

# sexual
tect_hutt_subset_enfa_spec_sexual <- enfa_hex_plot(tect_hutt_sexual_df_subset, marg = Mar, spec = Spe1, repro_mode = "Sexual")

ggsave("tect_hutt_subset_enfa_spec_sexual.png",
       plot = tect_hutt_subset_enfa_spec_sexual,
       device = "png",
       path = plot_path,
       dpi = "retina")

### 2) ENFA histogram
# asexual
hist(tect_hutt_asexual_enfa_subset)
title(main = "Asexual", adj = 0.7, line = -12)

png(filename = file.path(plot_path, "tect_hutt_subset_asexual_enfa_hist.png"))
hist(tect_hutt_asexual_enfa_subset)
title(main = "Asexual", adj = 0.7, line = -12)
dev.off()


# sexual
hist(tect_hutt_sexual_enfa_subset)
title(main = "Sexual", adj = 0.7, line = -12)

png(filename = file.path(plot_path, "tect_hutt_subset_sexual_enfa_hist.png"))
hist(tect_hutt_sexual_enfa_subset)
title(main = "Sexual", adj = 0.7, line = -12)
dev.off()

### 3) PCA of total e-space
tect_hutt_subset_total_pca <- total_climate_pca_plot(bg_env = tect_hutt_bg_env_subset, locs = tect_loc, genus = "Tectarchus", species = "huttoni")

tect_hutt_subset_total_pca

ggsave("tect_hutt_subset_total_pca.png",
       plot = tect_hutt_subset_total_pca,
       device = "png",
       path = plot_path,
       dpi = "retina")

# save enfa objects and remove them from the environment. They take up a lot of memory.
saveRDS(tect_hutt_sexual_enfa, file = file.path(obj_path, "tect_hutt_sexual_enfa.RDS"))
rm(tect_hutt_sexual_enfa)

saveRDS(tect_hutt_asexual_enfa, file = file.path(obj_path, "tect_hutt_asexual_enfa.RDS"))
rm(tect_hutt_asexual_enfa)

saveRDS(tect_hutt_sexual_enfa_subset,
        file = file.path(obj_path, "tect_hutt_subset_sexual_enfa.RDS"))
rm(tect_hutt_sexual_enfa_subset)

saveRDS(tect_hutt_asexual_enfa_subset,
        file = file.path(obj_path, "tect_hutt_subset_asexual_enfa.RDS"))
rm(tect_hutt_asexual_enfa_subset)

Stabilityy

We’re also interested in seeing if asexual populations live in more unstable climatic areas relative to sexual populations.

tect_hutt_locs_asexual <- tect_hutt_loc %>% 
  mutate(lat_long = str_c(latitude, longitude, sep = "_")) %>% 
  filter(reproductive_mode == "asexual",
         species == "huttoni",
         !duplicated(lat_long)) %>% 
  dplyr::select(longitude, latitude)

tect_hutt_locs_sexual <- tect_hutt_loc %>% 
  mutate(lat_long = str_c(latitude, longitude, sep = "_")) %>% 
  filter(reproductive_mode == "sexual",
         species == "huttoni",
         !duplicated(lat_long)) %>% 
  dplyr::select(longitude, latitude)

precip_asexual_tect_hutt <- raster::extract(precip_stability_scaled, tect_hutt_locs_asexual) %>% 
  enframe(name = NULL, value = "precip_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

precip_sexual_tect_hutt <- raster::extract(precip_stability_scaled, tect_hutt_locs_sexual) %>% 
  enframe(name = NULL, value = "precip_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

precip_df_tect_hutt <- bind_rows(precip_asexual_tect_hutt, precip_sexual_tect_hutt)

readr::write_csv(precip_df_tect_hutt, 
                 file.path(spread_path, "tect_hutt_precip_stability_df.csv"))

tect_hutt_precip_stability_plot <- ggplot(data = precip_df_tect_hutt, aes(x = reproductive_mode, y = precip_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

tect_hutt_precip_stability_plot

ggsave("tect_hutt_precip_stability.png",
       plot = tect_hutt_precip_stability_plot,
       device = "png",
       path = plot_path,
       dpi = "retina")

temp_asexual_tect_hutt <- raster::extract(temp_stability_scaled, tect_hutt_locs_asexual) %>% 
  enframe(name = NULL, value = "temp_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

temp_sexual_tect_hutt <- raster::extract(temp_stability_scaled, tect_hutt_locs_sexual) %>% 
  enframe(name = NULL, value = "temp_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

temp_df_tect_hutt <- bind_rows(temp_asexual_tect_hutt, temp_sexual_tect_hutt)

readr::write_csv(temp_df_tect_hutt, 
                 file.path(spread_path, "tect_hutt_temp_stability_df.csv"))

tect_hutt_temp_stability_plot <- ggplot(data = temp_df_tect_hutt, aes(x = reproductive_mode, y = temp_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

tect_hutt_temp_stability_plot

ggsave("tect_hutt_precip_stability.png",
       plot = tect_hutt_temp_stability_plot,
       device = "png",
       path = plot_path,
       dpi = "retina")

overall_asexual_tect_hutt <- raster::extract(overall_stability_scaled, tect_hutt_locs_asexual) %>% 
  enframe(name = NULL, value = "overall_stability_scaled") %>% 
  mutate(reproductive_mode = "asexual")

overall_sexual_tect_hutt <- raster::extract(overall_stability_scaled, tect_hutt_locs_sexual) %>% 
  enframe(name = NULL, value = "overall_stability_scaled") %>% 
  mutate(reproductive_mode = "sexual")

overall_df_tect_hutt <- bind_rows(overall_asexual_tect_hutt, overall_sexual_tect_hutt)

readr::write_csv(overall_df_tect_hutt, 
                 file.path(spread_path, "tect_hutt_overall_stability_df.csv"))

tect_hutt_overall_stability_plot <- ggplot(data = overall_df_tect_hutt, aes(x = reproductive_mode, y = overall_stability_scaled, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

tect_hutt_overall_stability_plot

ggsave("tect_hutt_overall_stability.png",
       plot = tect_hutt_overall_stability_plot,
       device = "png",
       path = plot_path,
       dpi = "retina")

Seasonality

seas_precip_asexual_tect_hutt <- raster::extract(w$bio15, tect_hutt_locs_asexual) %>% 
  enframe(name = NULL, value = "precip_seasonality") %>% 
  mutate(reproductive_mode = "asexual")

seas_precip_sexual_tect_hutt <- raster::extract(w$bio15, tect_hutt_locs_sexual) %>% 
  enframe(name = NULL, value = "precip_seasonality") %>% 
  mutate(reproductive_mode = "sexual")

seas_precip_df_tect_hutt <- bind_rows(seas_precip_asexual_tect_hutt, seas_precip_sexual_tect_hutt)

readr::write_csv(seas_precip_df_tect_hutt, 
                 file.path(spread_path, "tect_hutt_precip_seasonality_df.csv"))

tect_hutt_precip_seasonality_plot <- ggplot(data = seas_precip_df_tect_hutt, aes(x = reproductive_mode, y = precip_seasonality, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

tect_hutt_precip_seasonality_plot

ggsave(
  "tect_hutt_precip_seasonality.png",
  plot = tect_hutt_precip_seasonality_plot,
  device = "png",
  path = plot_path,
  dpi = "retina"
)

seas_temp_asexual_tect_hutt <- raster::extract(w$bio4, tect_hutt_locs_asexual) %>% 
  enframe(name = NULL, value = "temp_seasonality") %>% 
  mutate(reproductive_mode = "asexual")

seas_temp_sexual_tect_hutt <- raster::extract(w$bio4, tect_hutt_locs_sexual) %>% 
  enframe(name = NULL, value = "temp_seasonality") %>% 
  mutate(reproductive_mode = "sexual")

seas_temp_df_tect_hutt <- bind_rows(seas_temp_asexual_tect_hutt, seas_temp_sexual_tect_hutt)

readr::write_csv(seas_temp_df_tect_hutt, 
                 file.path(spread_path, "tect_hutt_temp_seasonality_df.csv"))

tect_hutt_temp_seasonality_plot <- ggplot(data = seas_temp_df_tect_hutt, aes(x = reproductive_mode, y = temp_seasonality, color = reproductive_mode)) +
  geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
  geom_jitter(width = 0.2, alpha = 0.6) +
  scale_color_viridis_d(option = "magma") +
  theme_dark()

tect_hutt_temp_seasonality_plot

ggsave(
  "tect_hutt_temp_seasonality.png",
  plot = tect_hutt_temp_seasonality_plot,
  device = "png",
  path = plot_path,
  dpi = "retina"
)

Move output

Put all output into species-specific subfolders.

Convenience scripts

These scripts aren’t crucial for reproducing this analysis, but were helpful for various reasons. Some of these have hard-coded paths and such, so no guarantees for use right out of the box.

This was a script I used to take full chelsa files, crop them to New Zealand extent, and write them to a file with my personal computer. I don’t have much memory, so unzipping to a temporary directory, then deleting the directory to free up space for the large files worked.

## get chelsa data
chelsa_folder <- "/Users/connorfrench/Dropbox/Old_Mac/climate-data/chelsa_30s_bio"
zip_files <- list.files(chelsa_folder, full.names = TRUE)

# using the Unarchiver commandline tools for Mac to unzip the 7zip chelsa layers. Regular unzip() does not work with 7z zipped files
for (file in zip_files) {
  # set temp directory
  tempd <- tempdir()
  system(paste("unar", file, "-o", tempd))
  r <- raster(list.files(tempd, pattern = "*.tif", full.names = TRUE)) %>%
    crop(extent(166, 179, -48, -34))
  writeRaster(r, filename = paste0("~/Desktop/", list.files(tempd, pattern = "*.asc")), format = "ascii")
  unlink(tempd, recursive = TRUE)
}